ADMB Documentation  11.1.2495
 All Classes Files Functions Variables Typedefs Friends Defines
multifan.cpp
Go to the documentation of this file.
00001 #include "statsLib.h"
00002 
00014 dvariable multifan(const dmatrix oprop,const dvar_matrix& pprop, const int& Nsamp)
00015 { //Vivian Haist.
00016     dvariable extra=0.1/14.;
00017     dvar_matrix resid=elem_div((oprop-pprop),sqrt((elem_prod(pprop,1.-pprop)+extra)/Nsamp));
00018     return sum(0.5*log(elem_prod(pprop,1. -pprop)+extra) -log(mfexp(-0.5*elem_prod(resid,resid))+0.01));
00019 }
00020 
00021 dvariable multifan(const int& n, const dmatrix obs, const dvar_matrix pred,double& nef)
00022 {
00023   int A=obs.colmax()-obs.colmin()+1;
00024   //dvar_matrix xi=(elem_prod(1.-pred,pred)+0.1/A)/n; //variance from Fourniers paper.
00025   dvar_matrix xi=(elem_prod(1.-obs,obs)+0.1/A)/n;  //variance from the multifanCL manual.
00026   dvar_matrix resid=obs-pred;
00027   nef=value(sum(elem_prod(1.-pred,pred))/sum(elem_prod(resid,resid)));
00028   return sum(0.5*log(2.*M_PI*xi)-log(mfexp(-0.5*elem_div(elem_prod(resid,resid),xi))+0.01));
00029 }
00030 
00031 dvariable multifan(const double& s,const dvector obsQ,const dvar_vector& preQ, double& nmle)
00032 {
00033   //using Fournier's robust likelihood for length frequency data.
00034   //s is the sample size
00035   //neff is the sample size limit  This seems to be fucked...
00036   //RETURN_ARRAYS_INCREMENT();
00037   dvariable like;
00038   dvariable tau;
00039   int lb=obsQ.indexmin();
00040   int nb=obsQ.indexmax();
00041 
00042   dvar_vector epsilon(lb,nb);
00043   dvar_vector Q=obsQ/sum(obsQ);
00044   dvar_vector Qhat=preQ/sum(preQ);
00045 
00046   //dvariable nmle;   //effective sample size
00047   nmle=value(sum(elem_prod(Qhat,1.-Qhat))/norm2(Q-Qhat));
00048   cout<<nmle<<endl;
00049   tau=1./s;
00050   epsilon=elem_prod(1.-Qhat,Qhat);
00051 
00052   like=0.5*sum(log(2.*M_PI*(epsilon+0.1/nb)))+nb*log(sqrt(tau));
00053   like+= -1.*sum(log(mfexp(-1.*elem_div(square(Q-Qhat),2.*tau*(epsilon+0.1/nb)))+0.01));
00054   //RETURN_ARRAYS_DECREMENT();
00055   return like;
00056 }
00057