Revision 401

trunk/examples/admb-re/priors_test/priors_user_guide.txt (revision 401)
1
ADMB new feature: Prior and likelihood section
2

  
3

  
4

  
5
The priors and likelihood sections was aim at help converting the winbugs model to admb. So the basic format is the same as winbugs  x ~ dnorm(mu,sigma); You can have any spaces between the ~ and within the argument lists. Here the dnorm is just any function name which taken on -log scale of the distribution, admb will provide some, you can has your own function or use the contribute library in contrib folder. One thing need notice is the first argument of actual function will be put on the left side of ~, all the remaining will be on right side of ~, such as your real function definition dnorm(x, mu,sigma), when you use ~, the first argument x will be on the left side of ~, the right side of ~ will be dnorm(mu,sigma), it doesn?t need to be x, but should has the same data type as x.  
6

  
7
The section definition can be PRIOR_SECTION or PRIORS_SECTION, LIKELIHOOD_SECTION or LIKELIHOODS_SECTION. The user can choose to use any of these section, which allowed to be put before procedure_section, if you choose to use both of them, the basic order will be priors section, likelihood section, and procedure section. For random effect model, the preliminary_calcs section only be allowed after likelihood section and before procedure section. Also for random effect, there are other limits too, one is the initialization section should be before parameter section if you use it, the runtime section only can be put after either the report section or preliminary section, I suggest you put after the report section.
8

  
9
The priors, likelihood sections allow to use any valid c++ code inside, it is like procedure_section, so don't forget the ; at the end of line except the if,else, for, while etc, and make sure they are valid c++ code. We will constrain the use of prior section only for parameters which are prefix init_ from parameter section. For the likelihood section, we constrain the use for constant data type, which from data_section only. If not, it throws the warning. But it is just remind warning, will not affect your final results. We also initialize the objective function value as 0 at the beginning of the procedure section, so you don't need to reinitialize, although it is not harm if you still want to do so. Remember every ~ you defined will be count into the objective function automatically, so use with some caution.
10

  
11
We also defined two variables, prior_function_value and likelihood_function_value, each sum up the negative log likelihood function value for their respective section, and will be added to objective function value automatically. Actually these two sections (including any C++ code there) will be appended at the end of procedure section.  You can get the value by referring the above name in your report section, such as report<<prior_function_value<<endl; 
12

  
13
To avoid the layout trouble of different sections, I would suggest you follow this order for random effect model, data section, initialization, parameter, priors, likelihood, preliminary, procedure, [separable function, function section,] between phase section, report section, [runtime section, final section, top of main, global section].  There is no order for sections within [ ], and the last two also can be put on the very top before data section.
14

  
15
For regular model, there is not much limit on the order, basically you follow data, parameter, priors, likelihood, procedure, report section, you can put initialization around the parameter, preliminary around the priors or likelihood, all other sections can be put almost anywhere.
16

  
trunk/examples/admb-re/priors_test/ratre.tpl (revision 401)
1
GLOBALS_SECTION
2
  #include "priors.cpp" 
3
  // append to mcmc output file
4
  long mcmc_iter = 0; 
5
  adstring mcFileName="admb1.mc.dat"; //put your mcmc saving file name
6
  
7
  //*** if your admb not build from the latest trunk, you need uncomment out following functions ***
8
  //df1b2variable operator -(double x, const df1b2variable& y) {return x+(-1.0*y);}
9
  //df1b2variable operator -(const df1b2variable& x,double y) {return (-1.0*y)+x;}
10

  
11

  
12
TOP_OF_MAIN_SECTION
13
  arrmblsize = 90000000L;
14
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(3000000);
15
  gradient_structure::set_CMPDIF_BUFFER_SIZE(2000000);
16
  gradient_structure::set_MAX_NVAR_OFFSET(8000);
17

  
18

  
19
DATA_SECTION 
20
  !!ad_comm::change_datafile_name("rat.dat");
21
  init_int N 
22
  init_int T 
23
  init_number xbar 
24
  init_vector x(1,T) 
25
  init_matrix y(1,N,1,T)
26
  //!! cout<<" ** x"<<endl<<x<<endl<<endl<<" ** y"<<endl<<y<<endl;// exit(8);
27

  
28
INITIALIZATION_SECTION
29
  tauc  1 
30
  alphac  150
31
  alphatau  1
32
  betac  10
33
  betatau 1
34

  
35
PARAMETER_SECTION 
36
  init_bounded_number tauc(0.00001,100,3) 
37
  init_number alphac    
38
  init_bounded_number alphatau(0.00001,100,2) 
39
  init_number betac 
40
  init_bounded_number betatau(0.00001,100,2) 
41
  random_effects_vector alpha(1,N) 
42
  random_effects_vector beta(1,N) 
43
  matrix mu(1,N,1,T) 
44
  sdreport_number alpha0   //user interest
45
  objective_function_value f 
46

  
47
PRIORS_SECTION
48
  tauc     ~ nllGamma ( 0.001 , 0.001); 
49
  alphac   ~ nllNormal2 ( 0.0 ,  0.000001);  
50
  alphatau ~ nllGamma (0.001 , 0.001); 
51
  betac    ~ nllNormal2( 0.0,   0.000001); 
52
  betatau  ~ nllGamma(0.001, 0.001); 
53
  //cout<<"1f ="<<prior_function_value<<"  "<< likelihood_function_value<<"  "<<f<<endl;
54

  
55
LIKELIHOOD_SECTION
56
  for(int i=1;i<=N;i++){ 
57
    alpha(i) ~ nllNormal2(alphac,alphatau);
58
    beta(i) ~ nllNormal2(betac,betatau);
59
    for(int j=1;j<=T;j++) 
60
    { 
61
      y(i,j) ~ nllNormal2(mu(i,j), tauc);
62
    }
63
  }
64
  //cout<<"2f ="<<prior_function_value<<"  "<< likelihood_function_value<<"  "<<f<<endl;
65
  //exit(8);
66

  
67
PRELIMINARY_CALCS_SECTION
68
  //cout<<T<<endl;
69

  
70
PROCEDURE_SECTION 
71
  //cout<<"alphac, betac = "<<alphac<<" "<<betac<<endl; //exit(8);
72
  for(int i=1;i<=N;i++) 
73
  { 
74
    for(int j=1;j<=T;j++) 
75
    { 
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff