// Extended version of: Orange data used in Ch. 8.2, Pinheiro & Bates (2000) // Correlated random effects 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) // Group indicator init_int m // Number of parameters in nonlinear regression model vector ofset(1,3); PARAMETER_SECTION init_bounded_vector beta(1,3,-50,200,1) // Fixed effects parameters init_bounded_number log_sigma(-6.0,5.0,2) // log(residual variance) init_bounded_vector log_sigma_u(1,3,-5.,5.,2) // 0.5*log(variance component) init_bounded_vector off_diag(1,3,-5.,5.,-1) // Off-diagonal elements in Cholesky factor random_effects_matrix u(1,M,1,3,2) objective_function_value g PRELIMINARY_CALCS_SECTION cout << setprecision(4); ofset(1) = 192.0; ofset(2) = 726.0; ofset(3) = 356.0; GLOBALS_SECTION #include //#include PROCEDURE_SECTION int i,ii; g = 0.0; ii = 0; for(i=1;i<=(int) M;i++) { fit_individual_tree(beta,u(i),i,ii,log_sigma,log_sigma_u,off_diag); } SEPARABLE_FUNCTION void fit_individual_tree(const dvar_vector& beta, const dvar_vector& u, const int& i, int& ii, const dvariable& log_sigma, const dvar_vector& log_sigma_u, const dvar_vector& off_diag) int j; dvar_matrix L(1,3,1,3); #ifdef OPT_LIB L.initialize(); #endif L(1,1) = 1.0; int k = 1; for(int i1=2;i1<=3;i1++) { L(i1,i1) = 1.0; for(j=1;j<=i1-1;j++) L(i1,j) = off_diag(k++); L(i1)(1,i1) /= norm(L(i1)(1,i1)); } //dvector offset("{192.0,726.0,356.0}"); // This does not appear to work. Why? dvar_vector a(1,3); //a = beta + ofset + elem_prod(mfexp(log_sigma_u),L*u); a = beta + elem_prod(mfexp(log_sigma_u),L*u); a(1) += 192.0; a(2) += 726.0; a(3) += 356.0; // Random effects contribution g -= -3*0.91893853 - .5*norm2(u); dvariable tmp, f; for(j=1;j<=ngroup(i);j++) { ii++; f = a(1)/(1+mfexp(-(t(ii)-a(2))/a(3))); tmp = y(ii) - f; tmp /= mfexp(log_sigma); g -= -0.91893853 - log_sigma - 0.5*square(tmp); } REPORT_SECTION report << beta + ofset << endl; report << exp(log_sigma) << endl; report << exp(log_sigma_u) << 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);