-
Notifications
You must be signed in to change notification settings - Fork 1
/
FG_MLR.stan
66 lines (62 loc) · 1.19 KB
/
FG_MLR.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
functions {
real FG_lpdf(vector y, real w1, real scale1, real scale2, int N) {
vector[N] z1;
vector[N] z2;
real p1;
vector[N] p2;
real p3;
vector[N] p4;
real output;
real w2;
w2 = 1.0 - w1;
z1 = y/scale1;
z2 = y/scale2;
p1 = w1/scale1;
p2 = exp(-z1-exp(-z1));
p3 = w2/scale2;
p4 = exp(z2-exp(z2));
output = sum(log(p1*p2 + p3*p4));
return output;
}
}
data {
int<lower=0> N;
int<lower=0> P;
vector[N] y;
matrix[N,P] X;
}
parameters {
real<lower=0, upper=1> w1;
real<lower=0> scale1;
real<lower=0> scale2;
vector[P] beta;
real alpha;
}
model {
alpha ~ normal(0,10^2);
beta ~ normal(0,10^2);
scale1 ~ inv_gamma(1,1);
scale2 ~ inv_gamma(1,1);
target += FG_lpdf(y-alpha-X*beta | w1,
scale1,
scale2,
N);
}
generated quantities {
vector[N] log_lik;
vector[N] ystar;
int z;
real x1;
real x2;
for (i in 1:N) {
z = bernoulli_rng(w1);
x1 = gumbel_rng(0,scale1);
x2 = -gumbel_rng(0,scale2);
ystar[i] = x1*z + x2*(1-z);
ystar[i] = ystar[i] + alpha + X[i,]*beta;
log_lik[i] = FG_lpdf(rep_vector(y[i]-alpha-X[i,]*beta,1) | w1,
scale1,
scale2,
1);
}
}