Mercurial > forge

changeset 11345:e75a08c7bfce octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
increased threshold for using series form of log likelihood function in gevlike
author nir-krakauer
date Wed, 02 Jan 2013 22:04:17 +0000
parents 4ebe6ca89b4f
children 558eb88d81b1
files main/statistics/inst/gevlike.m
diffstat 1 files changed, 12 insertions(+), 15 deletions(-) [+]
line wrap: on
line diff
--- a/main/statistics/inst/gevlike.m	Wed Jan 02 18:02:36 2013 +0000
+++ b/main/statistics/inst/gevlike.m	Wed Jan 02 22:04:17 2013 +0000
@@ -33,7 +33,7 @@
 ## @item
 ## @var{nlogL} is the negative log-likelihood.
 ## @item
-## @var{Grad} is the 3 by 1 gradient vector (first derivative of the log likelihood with respect to the parameter values)
+## @var{Grad} is the 3 by 1 gradient vector (first derivative of the negative log likelihood with respect to the parameter values)
 ## @item
 ## @var{ACOV} is the 3 by 3 Fisher information matrix (second derivative of the negative log likelihood with respect to the parameter values)
 ## 
@@ -76,18 +76,15 @@
   mu = params(3);
   
   #calculate negative log likelihood
-  if isargout(1)
-    [nll, k_terms] = gevnll (data, k, sigma, mu);
-    nlogL = sum(nll(:));
-  endif
+  [nll, k_terms] = gevnll (data, k, sigma, mu);
+  nlogL = sum(nll(:));
   
   #optionally calculate the first and second derivatives of the negative log likelihood with respect to parameters
-  if isargout(2)
-  	 [Grad, kk_terms] = gevgrad (data, k, sigma, mu, k_terms);
-  endif
-    
-  if isargout(3)
-  	ACOV = gevfim (data, k, sigma, mu, k_terms, kk_terms);
+  if nargout > 1
+  	 [Grad, kk_terms] = gevgrad (data, k, sigma, mu, k_terms);    
+    if nargout > 2
+    	ACOV = gevfim (data, k, sigma, mu, k_terms, kk_terms);
+    endif
   endif
 
 endfunction
@@ -104,11 +101,11 @@
     nlogL = exp(-a) + a + log(sigma);
   else
     aa = k .* a;
-    if min(abs(aa)) < 1E-4 && max(abs(aa)) < 0.5 #use a series expansion to find the log likelihood more accurately when k is small
+    if min(abs(aa)) < 1E-3 && max(abs(aa)) < 0.5 #use a series expansion to find the log likelihood more accurately when k is small
       k_terms = 1; sgn = 1; i = 0;
       while 1
         sgn = -sgn; i++;
-        newterm = sgn * (aa .^ i) / (i + 1);
+        newterm = (sgn  / (i + 1)) * (aa .^ i);
         k_terms = k_terms + newterm;
         if max(abs(newterm)) <= eps
           break
@@ -168,7 +165,7 @@
   kk_terms = 0.5; sgn = 1; i = 0;
   while 1
     sgn = -sgn; i++;
-    newterm = sgn * (aa .^ i) * ((i + 1) / (i + 2));
+    newterm = (sgn * (i + 1) / (i + 2)) * (aa .^ i);
     kk_terms = kk_terms + newterm;
     if max(abs(newterm)) <= eps
       break
@@ -257,7 +254,7 @@
   kkk_terms = 2/3; sgn = 1; i = 0;
   while 1
     sgn = -sgn; i++;
-    newterm = sgn * (aa .^ i) * ((i + 1) * (i + 2) / (i + 3));
+    newterm = (sgn * (i + 1) * (i + 2) / (i + 3)) * (aa .^ i);
     kkk_terms = kkk_terms + newterm;
     if max(abs(newterm)) <= eps
       break