forked from leaorsini/pseudo_gmm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Bayesian_GMM_on_simulated_pseudo_ar1_wcm_vector.stan
129 lines (101 loc) · 3.62 KB
/
Bayesian_GMM_on_simulated_pseudo_ar1_wcm_vector.stan
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
data {
int <lower=0> n; // nb of individuals
int <lower=0> K; // nb of clusters
real N; // nb of individuals
int< lower = 0> np; // np of parameters
matrix[n*K, np] X; // covariable matrix variable
vector[n*K] Y; // outcome variable
}
parameters {
vector[np] beta; // vector of parameters
}
model {
int loc1 = 0;
int loc2 = 0;
vector [K] yi;
matrix [K, np] xi; // covariable matrix variable
vector [K] ui;
matrix[K,K] fui_dev;
matrix[K,K] M1;
matrix[K,K] M2;
matrix[np,K] Di;
vector[2*np] gi;
vector[2*np] U = rep_vector(0,2*np);
matrix[2*np,2*np] C = rep_matrix(0, 2*np, 2*np);
matrix[2*np,2*np] Sigma;
//priors:
beta~ normal(0, sqrt(10)); //normal(0, 1);
M1[1,1] = 0.0; M1[1,2] = 1.0; M1[1,3] = 0.0; M1[1,4] = 0.0; M1[1,5] = 0.0;
M1[2,1] = 1.0; M1[2,2] = 0.0; M1[2,3] = 1.0; M1[2,4] = 0.0; M1[2,5] = 0.0;
M1[3,1] = 0.0; M1[3,2] = 1.0; M1[3,3] = 0.0; M1[3,4] = 1.0; M1[3,5] = 0.0;
M1[4,1] = 0.0; M1[4,2] = 0.0; M1[4,3] = 1.0; M1[4,4] = 0.0; M1[4,5] = 1.0;
M1[5,1] = 0.0; M1[5,2] = 0.0; M1[5,3] = 0.0; M1[5,4] = 1.0; M1[5,5] = 0.0;
//M2[1,1] = 1.0; M2[1,2] = 0.0; M2[1,3] = 0.0; M2[1,4] = 0.0; M2[1,5] = 0.0;
//M2[2,1] = 0.0; M2[2,2] = 0.0; M2[2,3] = 0.0; M2[2,4] = 0.0; M2[2,5] = 0.0;
//M2[3,1] = 0.0; M2[3,2] = 0.0; M2[3,3] = 0.0; M2[3,4] = 0.0; M2[3,5] = 0.0;
//M2[4,1] = 0.0; M2[4,2] = 0.0; M2[4,3] = 0.0; M2[4,4] = 0.0; M2[4,5] = 0.0;
//M2[5,1] = 0.0; M2[5,2] = 0.0; M2[5,3] = 0.0; M2[5,4] = 0.0; M2[5,5] = 1.0;
for(i in 1:n){
// set start location for next xi
loc1 = loc2+1;
loc2 = loc1+K-1;
yi = Y[loc1:loc2];
xi = X[loc1:loc2,];
ui = exp(-exp(xi*beta));
fui_dev = diag_matrix(-exp(xi*beta))*diag_matrix(ui);
Di = xi'*fui_dev;
gi[1:np] = (1/N)*Di*(yi-ui);
gi[(np+1):2*np] = (1/N)*Di*M1*(yi-ui);
//gi[(2*np+1):3*np] = (1/N)*Di*M2*(yi-ui);
C = C + gi*gi';
U = U + gi;
}
Sigma = C - (1/N)*U*U';
//pseudolikelihood for GMM :
target += -0.5*U'/Sigma*U;
}
generated quantities {
real loglik;
{
int loc1 = 0;
int loc2 = 0;
vector [K] yi;
matrix [K, np] xi; // covariable matrix variable
vector [K] ui;
matrix[K,K] fui_dev;
matrix[np,K] Di;
matrix[K,K] M1;
matrix[K,K] M2;
vector[2*np] gi;
vector[2*np] U = rep_vector(0,2*np);
matrix[2*np,2*np] C = rep_matrix(0, 2*np, 2*np);
matrix[2*np,2*np] Sigma;
M1[1,1] = 0.0; M1[1,2] = 1.0; M1[1,3] = 0.0; M1[1,4] = 0.0; M1[1,5] = 0.0;
M1[2,1] = 1.0; M1[2,2] = 0.0; M1[2,3] = 1.0; M1[2,4] = 0.0; M1[2,5] = 0.0;
M1[3,1] = 0.0; M1[3,2] = 1.0; M1[3,3] = 0.0; M1[3,4] = 1.0; M1[3,5] = 0.0;
M1[4,1] = 0.0; M1[4,2] = 0.0; M1[4,3] = 1.0; M1[4,4] = 0.0; M1[4,5] = 1.0;
M1[5,1] = 0.0; M1[5,2] = 0.0; M1[5,3] = 0.0; M1[5,4] = 1.0; M1[5,5] = 0.0;
//M2[1,1] = 1.0; M2[1,2] = 0.0; M2[1,3] = 0.0; M2[1,4] = 0.0; M2[1,5] = 0.0;
//M2[2,1] = 0.0; M2[2,2] = 0.0; M2[2,3] = 0.0; M2[2,4] = 0.0; M2[2,5] = 0.0;
//M2[3,1] = 0.0; M2[3,2] = 0.0; M2[3,3] = 0.0; M2[3,4] = 0.0; M2[3,5] = 0.0;
//M2[4,1] = 0.0; M2[4,2] = 0.0; M2[4,3] = 0.0; M2[4,4] = 0.0; M2[4,5] = 0.0;
//M2[5,1] = 0.0; M2[5,2] = 0.0; M2[5,3] = 0.0; M2[5,4] = 0.0; M2[5,5] = 1.0;
for(i in 1:n){
// set start location for next xi
loc1 = loc2+1;
loc2 = loc1+K-1;
yi = Y[loc1:loc2];
xi = X[loc1:loc2,];
ui = exp(-exp(xi*beta));
fui_dev = diag_matrix(-exp(xi*beta))*diag_matrix(ui);
Di = xi'*fui_dev;
gi[1:np] = (1/N)*Di*(yi-ui);
gi[(np+1):2*np] = (1/N)*Di*M1*(yi-ui);
//gi[(2*np+1):3*np] = (1/N)*Di*M2*(yi-ui);
C = C + gi*gi';
U = U + gi;
}
Sigma = C - (1/N)*U*U';
loglik = -0.5*U'/Sigma*U; // added to the output
}
}