[ADMB Users] adromb() ; numerical integration for RE models

[dave fournier]

On 12-11-13 09:16 AM, H. Skaug wrote:

I alkalized this approach would never work easily for the df1b2 stuff.
Better to put the necessary versions of the adromb and trapzd functions into
the  tpl file.  Need to get rid if the

    static s;

  line.  After tpl2rem need to get rid of extra
initialization lines near the bottom of the file.  This is a silly 
example derived from the
forestry example, but it should serve to show what needs to be done to
have Romberg integration in random effects files.  We really should have
a funnel_df1b2variable construction.



> Hi,
>
> Last week it was pointed out that adromb() does not work for models
> with random effects. Dave has kindly written code to fixed this (see
> attahced reredo.zip) but there must be something wrong
> in my settup, because it does not run properly for me (but it does
> for Dave), hence I do not want to submit to the code to the trunk.
>
> Dave is doing a big service to the community by providing code
> like this very rapidly. I thus feel sorry that I am not able
> to submit it. If somebody else are able to do it would be great.
>
> Hans

-------------- next part --------------
// Copyright (c) 2008, 2009, 2010 Regents of the University of California.
//
// ADModelbuilder and associated libraries and documentations are
// provided under the general terms of the "BSD" license.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2.  Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3.  Neither the name of the  University of California, Otter Research,
// nor the ADMB Foundation nor the names of its contributors may be used
// to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

DATA_SECTION
 int __i
 init_int nsteps 
 init_int k
 init_vector a(1,k+1)
 init_vector freq(1,k)
 int a_index;
 number sum_freq
!! sum_freq=sum(freq);
PARAMETER_SECTION
  init_bounded_number log_tau(-14,15,2)
  init_bounded_number log_nu(-15,4)
  init_bounded_number beta(.1,1.0,-1)
  init_bounded_number log_sigma(-5,3)
  likeprof_number tau
  sdreport_number nu
  sdreport_number sigma
  vector S(1,k+1)
  objective_function_value f
  random_effects_vector u(1,k+1)
INITIALIZATION_SECTION
  log_tau 0  
  beta 0.6666667 
  log_nu 0  
  log_sigma -2
PROCEDURE_SECTION
  tau=exp(log_tau);
  nu=exp(log_nu);
  sigma=exp(log_sigma);
   funnel_dvariable Integral;
   int i;
   for (i=1;i<=k+1;i++)
   {
     a_index=i;
     Integral=adromb(-3.0,3.0,nsteps);
     S(i)=Integral;
     //S(i)=h(0);
   }
   f=0.0;
   for (i=1;i<=k;i++)
   {
     dvariable ff=0.0;
     dvariable diff=posfun((S(i)-S(i+1))/S(i),.000001,ff);
     f-=freq(i)*log(1.e-50+S(i)*diff);
     f+=ff;
   }
   f+=sum_freq*log(1.e-50+S(1));
   f+=0.5*norm2(u);

FUNCTION   dvariable adromb(double a,double b,int ns)
  const double base = 4;
  const int JMAX=50;
  dvariable ss=0.0;
  int MAXN = min(JMAX, ns);
  dvar_vector s(1,MAXN+1);

  for(int j=1; j<=MAXN+1; j++)
  {
    s[j] = trapzd(a,b,j,ss);
  }

  for(int iter=1; iter<=MAXN+1; iter++)
  {
    for(int j=1; j<=MAXN+1-iter; j++)
    {
      s[j] = (pow(base,iter)*s[j+1]-s[j])/(pow(base,iter)-1);
    }
  }
  return s[1];

FUNCTION   dvariable trapzd(double a,double b,int n,const dvariable& s)
  double x,num_interval,hn;
  dvariable sum;
 // static dvariable s;
  static int interval;
  int j;
  model_parameters * ptr= (model_parameters *) mycast();
  if (n == 1) {
    interval=1;
    return (s=0.5*(b-a)*(h(a)+h(b)));
  } else {
    num_interval=interval;
    hn=(b-a)/num_interval;
    x=a+0.5*hn;
    for (sum=0.0,j=1;j<=interval;j++,x+=hn) sum += h(x);
    interval *= 2;
    s=0.5*(s+(b-a)*sum/num_interval);
    return s;
  }

FUNCTION dvariable h(const double z)
  dvariable tmp;
  tmp=exp(-.5*z*z + u(a_index)+ tau*(-1.+exp(-nu*pow(a(a_index),beta)*exp(sigma*z))) );  
  return tmp;

REPORT_SECTION
  report << "nsteps = " << std::setprecision(10) <<  nsteps << endl;
  report << "f = " << std::setprecision(10) <<  f << endl;
  report << "a" << endl << a << endl;
  report << "freq" << endl << freq << endl;
  report << "S" << endl << S << endl;
  report << "S/S(1)" << endl << std::ios::fixed << std::setprecision(6) << S/S(1) << endl;
  report << "tau "  << tau << endl; 
  report << "nu "  << nu << endl; 
  report << "beta "  << beta << endl; 
  report << "sigma "  << sigma << endl; 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: forest.dat
Type: application/x-ns-proxy-autoconfig
Size: 260 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121113/284f0d79/attachment.dat>