ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
dmvlogistic.cpp
Go to the documentation of this file.
00001 #include "statsLib.h"
00002 
00015 dvariable dmvlogistic(const dmatrix o, const dvar_matrix& p,dvar_matrix& nu, double& tau2,const double minp)
00016 { //returns the negative loglikelihood using the MLE for the variance
00017   /*
00018     This is a modified version of the dmvlogistic negative log likelihood
00019     where proportions at age less than minp are pooled into the consecutive
00020     age-classes.  See last paragraph in Appendix A of Richards, Schnute and
00021     Olsen 1997.
00022 
00023     NB minp must be greater than 0, otherwise this algorithm returns an
00024     error if one of the observed proportions is zero.
00025 
00026     -1) first count the number of observations for each year > minp
00027     -2) normalized observed and predicted age-proportions
00028     -3) loop over ages, and check if observed proportion is < minp
00029         -if yes, then add observed proprtion to bin k
00030         -if no then add observed proportion to bin k and increment
00031          bin k if k is currently less than the number of bins.
00032     -4) do the same grouping for the predicted proportions.
00033     -5) use ivector iiage to correctly assign residuals into nu
00034 
00035     FEB 8, 2011.  Fixed a bug in the variance calculation &
00036     likelihood scaling that was discovered at the 2011 Hake
00037     assessment STAR panel in Seattle.
00038   */
00039   RETURN_ARRAYS_INCREMENT();
00040   int i,j,k,n;
00041   int age_counts=0;
00042   int a = o.colmin();
00043   int A=o.colmax();
00044   double tiny=0.001/(A-a+1);
00045   int t=o.rowmin();
00046   int T=o.rowmax();
00047   dvariable tau_hat2;   //mle of the variance
00048   //dvar_matrix nu(t,T,a,A);
00049   nu.initialize();
00050 
00051   for(i=t; i<=T; i++)
00052   {
00053     n=0;
00054     dvector oo = o(i)/sum(o(i));
00055     dvar_vector pp = p(i)/sum(p(i));
00056 
00057     //count # of observations greater than minp (2% is a reasonable number)
00058     for(j=a;j<=A;j++)
00059       if(oo(j) > minp)n++;
00060 
00061     ivector iiage(1,n);
00062     dvector o1(1,n); o1.initialize();
00063     dvar_vector p1(1,n); p1.initialize();
00064     k=1;
00065     for(j=a;j<=A;j++)
00066     {
00067       if(oo(j)<=minp)
00068       {
00069         o1(k)+=oo(j);
00070         p1(k)+=pp(j);
00071       }
00072       else
00073       {
00074         o1(k)+=oo(j);
00075         p1(k)+=pp(j);
00076         if(k<=n)iiage(k)=j;   //ivector for the grouped residuals
00077         if(k<n) k++;
00078       }
00079     }
00080 
00081     //assign residuals to nu based on iiage index
00082     dvar_vector t1 = log(o1)-log(p1) - mean(log(o1)-log(p1));
00083 
00084     for(j=1;j<=n;j++)
00085       nu(i)(iiage(j))=t1(j);
00086 
00087     age_counts += n-1;
00088   }
00089   //Depricated  Wrong Variance & likelihood calculation.
00090   //tau_hat2 = 1./(age_counts*(T-t+1))*norm2(nu);
00091   //dvariable nloglike =(age_counts*(T-t+1))*log(tau_hat2);
00092 
00093   //Feb 8, 2011  Fixed variance & likelihood
00094   tau_hat2 = 1./(age_counts)*norm2(nu);
00095   dvariable nloglike =(age_counts)*log(tau_hat2);
00096   tau2=value(tau_hat2); //mle of the variance
00097   RETURN_ARRAYS_DECREMENT();
00098   return(nloglike);
00099 }
00100 
00101 //The following function has been depricated.
00102 dvariable dmvlogistic(const dmatrix o, const dvar_matrix& p, double& tau2)
00103 {
00104 
00105   //returns the negative loglikelihood using the MLE for the variance
00106   RETURN_ARRAYS_INCREMENT();
00107 
00108   int a = o.colmin();
00109   int A=o.colmax();
00110   double tiny=0.001/(A-a+1);
00111   int t=o.rowmin();
00112   int T=o.rowmax();
00113   dvariable tau_hat2;   //mle of the variance
00114   dvar_matrix nu(t,T,a,A);
00115   for(int i=t; i<=T; i++)
00116   {
00117     dvector o1=o(i)/sum(o(i));
00118     dvar_vector p1=p(i)/sum(p(i));
00119     //p1=log(p1)-mean(p1);
00120     //p1=log(p1)-mean(log(p1));
00121     //p1 = exp(p1)/sum(exp(p1));
00122     nu(i) = log(o1+tiny)-log(p1+tiny) - mean(log(o1+tiny)-log(p1+tiny));
00123   }
00124   tau_hat2 = 1./((A-a)*(T-t+1))*norm2(nu);
00125   dvariable nloglike =((A-a)*(T-t+1))*log(tau_hat2);
00126   tau2=value(tau_hat2); //mle of the variance
00127   RETURN_ARRAYS_DECREMENT();
00128   return(nloglike);
00129 }
00130