// Orange data used in Millar (2004) Aust NZ J. Stat Vol 46, p. 543-554. DATA_SECTION init_int n // Number of data points init_vector y(1,n) // Response vector init_vector t(1,n) // Primary covariate init_int M // Number of groups init_vector ngroup(1,M) // Number in each group init_int m // Number of parameters in nonlinear regression model PARAMETER_SECTION init_bounded_vector beta(1,3,-50,200,1) // Fixed effects parameters init_bounded_number log_sigma(-5.0,5.0,1) // log(residual variance) init_bounded_number log_sigma_u(-5,5,2) // 0.5*log(variance component) init_bounded_number log_sigma_v(-5,5,2) // 0.5*log(variance component) random_effects_vector u(1,M,2) // Unscaled random effects random_effects_vector v(1,7,2) objective_function_value g PRELIMINARY_CALCS_SECTION cout << setprecision(4); GLOBALS_SECTION #include PROCEDURE_SECTION int i,ii,j; dvar_vector a(1,3); // Basic model function parameters dvariable tmp, f; dvariable sigma = mfexp(log_sigma); a(1) = 192.0 + beta(1); a(2) = 726.0 + beta(2); a(3) = 356.0 + beta(3); g = 0.0; ii=0; for(i=1;i<=(int) M;i++) { for(j=1;j<=ngroup(i);j++) { ii++; g -= -log_sigma; f = (a(1)+u(i)+v(j))/(1+mfexp(-(t(ii)-a(2))/a(3))); tmp = y(ii) - f; tmp /= sigma; g -= -0.5*tmp*tmp; } } // Random effect contribution from u for(i=1;i<=M;i++) { g -= -(log_sigma_u); g -= -.5*(square(u(i)/mfexp(log_sigma_u))); } // Random effects contribution from v for(j=1;j<=7;j++) { g -= -(log_sigma_v); g -= -.5*(square(v(j)/mfexp(log_sigma_v))); } REPORT_SECTION //report << beta0+beta << endl; report << exp(log_sigma) << endl; report << exp(log_sigma_u) << endl; report << exp(log_sigma_v) << endl; TOP_OF_MAIN_SECTION arrmblsize = 40000000L; gradient_structure::set_GRADSTACK_BUFFER_SIZE(3000000); gradient_structure::set_CMPDIF_BUFFER_SIZE(200000); gradient_structure::set_MAX_NVAR_OFFSET(10000);