[ADMB Users] Fwd: lgam in 10.0

[dave fournier]

The third parameter (which I think of as tau) should be > 1.0
Something like this.


DATA_SECTION
         init_int nobs
         init_ivector ncones(1,nobs)
         init_vector DBH(1,nobs)
PARAMETER_SECTION
         init_bounded_number a(0.0,10.0)
         init_bounded_number b(0.0,10.0)
         init_bounded_number logk(-10.,10.0)
         sdreport_number k
         vector mu(1,nobs)
         objective_function_value f
PROCEDURE_SECTION
         k=1.0+exp(logk);
         f=0.0;
         mu=a*pow(DBH,b);
         for(int i=1; i<=nobs; i++)
         {
                 f-=log_negbinomial_density(ncones[i], mu[i], k);
         }
REPORT_SECTION
         double tmp=0.0;
         int i;
         for(i=1; i<=nobs; i++)
         {
            tmp+=value(square(ncones[i]-mu[i])/(k*mu[i]));
         }
         tmp/=nobs;
         report << tmp << endl;
         for(i=1; i<=nobs; i++)
         {
            report << ncones(i) << " " << mu(i) << "  "
<< square(ncones[i]-mu[i])/(k*mu[i])<< endl;
         }


Put in a little report to examine the fit (still raining but soon the 
hockey game.)

You have a few big outliers.
Maybe robustify the model as a mixture of NB's with the same means but a 
multiple of tau's (k;'s here).

        init_int nobs
         init_ivector ncones(1,nobs)
         init_vector DBH(1,nobs)
PARAMETER_SECTION
         init_bounded_number a(0.0,10.0)
         init_bounded_number b(0.0,10.0)
         init_bounded_number logk(-10.,10.0)
         init_bounded_number logpcon(-3.0,-0.5,3)
         init_bounded_number logkmult(-0.0,4.0,3)
         init_number kkludge(2)
         sdreport_number k
         sdreport_number kmult
         sdreport_number pcon
         vector mu(1,nobs)
         objective_function_value f
PROCEDURE_SECTION
         pcon=exp(logpcon);
         kmult=exp(logkmult);
         k=1.0+exp(logk);
         f=0.0;
         mu=a*pow(DBH,b);
         if (current_phase()<2)
         {
           for(int i=1; i<=nobs; i++)
           {
             f+=square(log( ncones(i)+1)-log(mu(i)+1.0));
           }
         }
         else
         {
           for(int i=1; i<=nobs; i++)
           {
             f-=log((1.0-pcon)*exp(log_negbinomial_density(ncones[i], 
mu[i], k))
               + pcon*exp(log_negbinomial_density(ncones[i], 
mu[i],kmult*k)));
           }
         }
         f+=square(kkludge);
REPORT_SECTION
         double tmp=0.0;
         int i;
         for(i=1; i<=nobs; i++)
         {
            tmp+=value(square(ncones[i]-mu[i])/(k*mu[i]));
         }
         tmp/=nobs;
         report << tmp << endl;
         for(i=1; i<=nobs; i++)
         {
            report << ncones(i) << " " << mu(i) << "  "
<< square(ncones[i]-mu[i])/(k*mu[i])<< endl;
         }



Fit this firt to get the means sort of right.

             f+=square(log( ncones(i)+1)-log(mu(i)+1.0));


Output reports par files are
# Number of parameters = 3  Objective function value = 1128.40  Maximum 
gradient component = 7.54442e-08
# a:
0.365011644781
# b:
2.24328259535
# logk:
3.277

and

# Number of parameters = 6  Objective function value = 1108.58  Maximum 
gradient component = 1.42155e-07
# a:
0.312960997627
# b:
2.32658393536
# logk:
2.70363429357
# logpcon:
-2.22117808691
# logkmult:
2.56571111968
# kkludge:
0.00000000000

So a big change in the LL for adding 3 parameters.

Hockey time ...