Chapter 3 A forestry model -- estimating the size distribution of wildfires

Section 1 Model description

This examples highlights two features of AD Model Builder, the use of a numerical integration routine within a statistical parameter estimation model and the use of the ad_begin_funnel mechanism to reduce the size of temporary file storage required. It also provides a performance comparison between AD Model Builder and Splus.

This problem investigates a model which predicts a relationship between the size and frequency of wildfires. It is assumed that the probability of observing a wildfire in size category is given by , where

is the number of widfires observed to lie in size category

the log-likelihood function for the problem is given by

where

is defined by the integral

The parameters , , , and are functions of the parameters of the original model, and don't have a simple interpretation. Fitting the model to data involves maximizing the above log-likelihood (\chapno.1). While the gradient can be calculated (in integral form), coding it is cumbersome. Numerically maximizing the log-likelihood without specifying the gradient is preferable.

The parameter is related to the fractal dimension of the perimeter of the fire. One hypothesis of interest is that which is related to hypotheses about the nature of the mechanism by which fires spread. The AD Model Builder code for the model follows.


DATA_SECTION
 int time0 
 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_number log_tau
  init_number log_nu
  init_number log_beta(2)
  init_number log_sigma
  sdreport_number tau
  sdreport_number nu
  sdreport_number sigma
  sdreport_number beta
  vector S(1,k+1)
  objective_function_value f
INITIALIZATION_SECTION
  log_tau 0  
  log_beta -.405465 
  log_nu 0  
  log_sigma -2
PROCEDURE_SECTION
  tau=exp(log_tau);
  nu=exp(log_nu);
  sigma=exp(log_sigma);
  beta=exp(log_beta);
   funnel_dvariable Integral;
   int i;
   for (i=1;i<=k+1;i++)
   {
     a_index=i;
     ad_begin_funnel();
     Integral=adromb(&model_parameters::h,-3.0,3.0,nsteps);
     S(i)=Integral;
   }
   f=0.0;
   for (i=1;i<=k;i++)
   {
     dvariable ff=0.0;
     // make the model stable for case when S(i)<=S(i+1)
     // we have to subrtract s(i+1) from S(i) first or roundoff will
     // do away with the 1.e-50.
     f-=freq(i)*log(1.e-50+(S(i)-S(i+1)));
     f+=ff;
   }
   f+=sum_freq*log(1.e-50+S(1));
FUNCTION dvariable h(const dvariable& z)
  dvariable tmp;
  tmp=exp(-.5*z*z + tau*(-1.+exp(-nu*pow(a(a_index),beta)*exp(sigma*z))) );  
  return tmp;
REPORT_SECTION
  int * pt=NULL;
  report << " elapsed time = "  << time(pt)-time0 << " seconds" << endl;
  report << "nsteps = " << setprecision(10) <<  nsteps << endl;
  report << "f = " << setprecision(10) <<  f << endl;
  report << "a" << endl << a << endl;
  report << "freq" << endl << freq << endl;
  report << "S" << endl << S << endl;
  report << "S/S(1)" << endl << setfixed << setprecision(6) << S/S(1) << endl;
  report << "tau "  << tau << endl; 
  report << "nu "  << nu << endl; 
  report << "beta "  << beta << endl; 
  report << "sigma "  << sigma << endl;

Section 2 The numerical integration routine

The statement


   Integral=adromb(&model_parameters::h,-3.0,3.0,nsteps);

invokes the numerical integration routine for the user-defined function h. The function must be defined in a FUNCTION subsection. It can have any name, must be defined to take a const dvariable& argument, and must return a dvariable. The values -3.0, 3.0 are the limits of integration (effectively

for this example). The integer argument nsteps determines how accurate the integration will be. Higher values of nsteps will be more accurate but greatly increase the amount of time necessary to fit the model. The basic strategy is to use a moderate value for nteps, such as 6, and then to increase this value to see if the parameter estimates change much.


FUNCTION dvariable h(const dvariable& z)

Section 3 Using the ad_begin_funnel routine to reduce the amount of temporary storage required

Numerical integration routines can be very computationally intensive, especially when they must be computed to great accuracy. Such computations will require a lot of temporary storage in AD Model Builder. Fortunately the output from such a routine is just one number, the value of the integral. In automatic differentiation terminology a long set of computations which produce just one number is known as a funnel. It is possbile to exploit the properties of such a funnel to greatly reduce the amount of temporary storage required. All that is necessary is to declare an object of type funnel_dvariable and to assign the results of the computation to it. At the beginning of the funnel a call to the function ad_begin_funnel is made. There is quite a bit of overhead associated with the funnel construction so it should not be used for very small calculations. However it is possible to put it in and test the program to see whether it runs more quickly or not. The following modifed code will produce exactly the same results, but without the funnel construction.


  dvariable Integral;   // change the definition of Integral
  int i;
  for (i=1;i<=k+1;i++)
  {
    a_index=i;
    // ad_begin_funnel();  // commment out this line
    Integral=adromb(&model_parameters::h,-3.0,3.0,nsteps);
    S(i)=Integral;
  }

If the funnel construction is used on a portion of code which is not a funnel, incorrect derivative values will be obtained. If this is suspected the funnel should be removed as in the above example and the model run again.

Section 4 Effect of the accuracy switch on the running time for numerical integration

The following report shows the amount of time required to run the model with a fxied value of

for different values of the parameter nsteps. For practical perposes a vlaue of nsteps=8 gives enough accuracy so that the model could be fit in about 6 seconds.


elapsed time = 2 seconds nsteps = 6 f = 629.9846518
tau 9.851110 nu 8.913479 beta 0.666667 sigma 1.885570

elapsed time = 2 seconds nsteps = 7 f = 629.9851092
tau 9.850213 nu 8.835066 beta 0.666667 sigma 1.882967

elapsed time = 6 seconds nsteps = 8 f = 629.9851223
tau 9.850227 nu 8.836769 beta 0.666667 sigma 1.883024

elapsed time = 6 seconds nsteps = 9 f = 629.9851222
tau 9.850226 nu 8.836769 beta 0.666667 sigma 1.883024

elapsed time = 14 seconds nsteps = 10 f = 629.9851222
tau 9.850226 nu 8.836769 beta 0.666667 sigma 1.883024

The corresponding times when beta was estimated in an extra phase of the minimization are given here. It as apparent tat the model parameters become unstable when beta is being estimated. Twice the log-likelihood difference is which is significant



elapsed time = 3 seconds nsteps = 6 f = 627.2919906
tau 20.729183 nu 427.816375 beta 0.180225 sigma 2.499445

elapsed time = 6 seconds nsteps = 7 f = 627.2952716
tau 21.868971 nu 80914.970724 beta 0.170392 sigma 4.232237

elapsed time = 17 seconds nsteps = 8 f = 627.297021
tau 22.858629 nu 2326271883.421848 beta 0.164749 sigma 7.653068

elapsed time = 62 seconds nsteps = 9 f = 627.2993787
tau 23.771061 nu 1652877622661391616.000000 beta 0.161073 sigma 14.451510

elapsed time = 123 seconds nsteps = 10 f = 627.3106333
tau 23.116097 nu 49753858778.636856 beta 0.159364 sigma 8.663666

elapsed time = 244 seconds nsteps = 11 f = 627.310624
tau 23.115275 nu 49009470510.133156 beta 0.159369 sigma 8.658643

Section 5 A comparison with Splus for the forestry model

The Splus minimizing routine nlminb was used to fit the model. Fitting the three parameter model with Splus required approximately 280 seconds compared to 6 seconds with AD Model Builder, so that AD Model Builder was approximately 45 times faster for this simple problem.

For the four parameter problem with beta estimated, the SPLUS routine exited after fourteen minutes and 30 seconds, reporting false convergence with a function value of 627.338.

The data for the example is


a
 0.04 0.1 0.2 0.4 0.8 1.6 3.2 6.4 12.8 25.6 51.2 102.4 204.8
freq
 167 84 61 29 19 17 4 4 1 0 1 1

where the first line contains the bounds for the size catagories and the second line contains the number of observations in each size category. The Splus code with fixed beta for the example is


obj.20<-
function(xvec)
{
#Objective for maxn in NLMINB   NB vector argument
 - llik.20(xvec[1], xvec[2], xvec[3])
}
llik.20<-
function(logtau, lognu, logsigma)
{
        tau<-exp(logtau)
        nu<-exp(lognu)
        sigma<-exp(logsigma)
print(tau)
print(nu)
print(sigma)
        llik <- 0
        for(i in 1:(length(freq)+1)) {
           Int[i]<-S.20(xa[i], tau, nu, sigma)
}      
        print(llik)
        for(i in 1:length(freq)) {
           llik <- llik + (freq[i] * (log(1.e-50+(Int[i]-Int[i+1])) 
               -log(1.e-50+Int[1])))
        }
        llik
}
S.20<-
function(da, tau, nu, sigma)
{
        results <- integrate(intgnd.20, -3, 3, TAU = tau, NU = nu, SIGMA =
                sigma, A = da)
        if(results$message != "normal termination")
                ans <- results$message
        else ans <- results$integral
        ans
}
intgnd.20<-
function(z, A, TAU, NU, SIGMA)
{
  exp(-0.5 * z^2 + TAU * (-1 +exp( - NU * A^(2/3) * exp(SIGMA * z))))
}

To run the example in Splus with the same initial values use the following values


  logtau 0  lognu 0  logsigma -2

The vector xa should contain the 13 a values while the vector freq should contain the 12 observed frequencies.

NEXT