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

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;

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)

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.

elapsed time = 2 seconds nsteps = 6 f = 629.9846518 tau 9.851110 nu 8.913479 beta 0.666667 sigma 1.885570elapsed 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

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 1where 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 -2The vector xa should contain the 13 a values while the vector freq should contain the 12 observed frequencies.