-
Notifications
You must be signed in to change notification settings - Fork 1
/
wt.tpl
179 lines (142 loc) · 4.1 KB
/
wt.tpl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
DATA_SECTION
init_number binit;
init_int syr;
init_int nyr;
init_int sage;
init_int nage;
init_matrix data(syr,nyr,sage,nage);
init_matrix obs_std(syr,nyr,sage,nage);
int byr;
!! byr = syr - (nage-sage+1);
!! cout<<byr<<endl;
!! cout<<syr<<endl;
!! cout<<nyr<<endl;
vector age(sage,nage);
!!age.fill_seqadd(sage,1);
INITIALIZATION_SECTION
log_winf 2.302585093;
log_vonb -1.049822124;
log_to -0.693147181;
b binit;
PARAMETER_SECTION
init_number log_winf(1);
init_number log_vonb(1);
init_number log_to(1);
init_bounded_number b(2,4,4);
init_bounded_number log_sd_year(-30,10,3);
init_bounded_number log_sd_cohort(-30,10,3);
// init_bounded_vector year_effect(syr,nyr,-5.0,5.0,2);
// init_bounded_vector cohort_effect(byr,nyr,-5.0,5.0,2);
random_effects_vector year_effect(syr,nyr,2);
random_effects_vector cohort_effect(byr,nyr,2);
objective_function_value f;
number winf;
number vonb;
number to;
number sd_year;
number sd_cohort;
number nll_year_effect;
number nll_cohort_effect;
vector wa(sage,nage);
matrix W(syr,nyr,sage,nage); // Predicted Weight
matrix E(syr,nyr,sage,nage); // Random effects
matrix R(syr,nyr,sage,nage); // Residuals
PROCEDURE_SECTION
winf = exp(log_winf);
vonb = exp(log_vonb);
to = - exp(log_to);
sd_year = mfexp(log_sd_year);
sd_cohort = mfexp(log_sd_cohort);
wa = winf * pow(1.0 - exp( -vonb*(age-to)),b );
//COUT(wa);
// Initial Cohort effects
int iyr = byr;
for( int j = sage; j <= nage; j++){
E(syr,j) = exp(sd_cohort * cohort_effect(iyr++) +
sd_year * year_effect(syr));
}
for( int i = syr+1; i <= nyr; i++){
for( int j = sage; j <= nage; j++){
if( j == sage ){
E(i,j) = exp(sd_cohort * cohort_effect(i));
//E(i,j)*= exp(sd_year * year_effect(i));
} else {
E(i,j) = E(i-1,j-1) * exp(sd_year * year_effect(i));
}
}
}
// predicted data where the weight increment for the vb model is given by:
// dw=%e^(-b*delta*k-a1*b*k) * ((%e^(delta*k+a1*k)-%e^(k*to))^b
// -%e^(b*delta*k)*(%e^(a1*k)-%e^(k*to))^b)*winf
double delta = 1; // annual time step.
dvariable t1,t2,t3,t4,t5,t6;
dvariable k,dw;
k = vonb;
W(syr) = elem_prod(wa,E(syr));
R(syr) = data(syr) - W(syr);
//cout<<W<<endl;
for( int i = syr+1; i <= nyr; i++){
for( int j = sage; j <= nage; j++){
if( j == sage ){
W(i,j) = wa(j) * E(i,j);
} else if(j > sage){
double a1 = double(j)-1.0;
t1 = exp(-b*delta*k-a1*b*k);
t2 = pow(exp(delta*k+a1*k)-exp(k*to),b);
t3 = exp(b*delta*k) * pow(exp(a1*k)-exp(k*to),b);
t5 = t2 - t3;
dw = winf * t1 * t5;
//cout << j << "\t" << dw<< endl;
W(i,j) = W(i-1,j-1) + dw * E(i,j);
}
// Residual
R(i,j) = data(i,j) - W(i,j);
}
}
// Compute the objective function
double log2pi = log(2.0*M_PI);
dvariable constant;
dvariable nloglike;
dvariable variance;
nloglike=0;
int n = size_count(age);
for( int i = syr+1; i <= nyr; i++){
for( int j = sage; j <= nage; j++){
variance = square(obs_std(i,j));
constant = 0.5 * log2pi + 0.5 * log(variance);
nloglike += constant + 1.0/(2*variance)*square(R(i,j));
}
}
// pdf for year effect
dvariable nll_year_effect = 0.0;
n = 0.5*size_count(year_effect);
nll_year_effect = n * log2pi + n * log(square(sd_year))
+ 0.5 * norm2(year_effect) / square(sd_year);
// nll_year_effect = 0.5 * norm2(year_effect);
// pdf for cohort effect
dvariable nll_cohort_effect = 0.0;
n = 0.5*size_count(cohort_effect);
nll_cohort_effect = n * log2pi + n * log(square(sd_cohort))
+ 0.5 * norm2(cohort_effect) / square(sd_cohort);
// nll_cohort_effect = 0.5 * norm2(cohort_effect);
f = nloglike + nll_year_effect + nll_cohort_effect;
//cout<<E<<endl;
//cout<<W<<endl;
//exit(1);
REPORT_SECTION
REPORT(winf);
REPORT(vonb);
REPORT(to);
REPORT(sd_year);
REPORT(sd_cohort);
REPORT(data);
REPORT(W);
REPORT(R);
REPORT(E);
GLOBALS_SECTION
#undef REPORT
#define REPORT(object) report << #object "\n" << setw(8) \
<< setprecision(6) << setfixed() << object << endl;
#undef COUT
#define COUT(object) cout << #object "\n" << setw(6) \
<< setprecision(5) << setfixed() << object << endl;