-
Notifications
You must be signed in to change notification settings - Fork 80
/
orange_big.cpp
57 lines (47 loc) · 1.08 KB
/
orange_big.cpp
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
// Scaled up version of the Orange Tree example (50,000 latent random variables)
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
DATA_INTEGER(n);
DATA_VECTOR(y);
DATA_VECTOR(t);
DATA_INTEGER(M);
DATA_FACTOR(ngroup);
DATA_INTEGER(multiply);
PARAMETER_VECTOR(beta);
PARAMETER(log_sigma);
PARAMETER(log_sigma_u);
PARAMETER_VECTOR(u);
Type sigma=exp(log_sigma);
Type sigma_u=exp(log_sigma_u);
ADREPORT(sigma);
ADREPORT(sigma_u);
using namespace density;
int i,j,k,ii;
parallel_accumulator<Type> g(this);
for(k=0;k< multiply;k++)
{
ii = 0;
for(i=0;i< M;i++)
{
// Random effects contribution
Type u1 = u[i+k*M];
g -= -.5*u1*u1;
vector<Type> a(3);
a[0] = 192.0 + beta[0] + sigma_u*u1;
a[1] = 726.0 + beta[1];
a[2] = 356.0 + beta[2];
Type tmp;
Type f;
for(j=0;j<ngroup(i);j++)
{
f = a[0]/(1+exp(-(t[ii]-a[1])/a[2]));
tmp = (y[ii] - f)/sigma;
g -= -log_sigma - 0.5*tmp*tmp;
ii++;
}
}
}
return g;
}