-
Notifications
You must be signed in to change notification settings - Fork 81
/
stan_mnl_altspecific.stan
137 lines (95 loc) · 3.51 KB
/
stan_mnl_altspecific.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
129
130
131
132
133
134
135
136
137
data {
int K; # number of choices
int N; # number of individuals
int D; # number of indiv specific vars
int G; # number of alt specific vars
int T; # number of alt constant vars
int y[N*K]; # choices
vector[N*K] choice; # choice made (logical)
matrix[N,D] X; # data for indiv specific effects
matrix[N*K, G] Y; # data for alt specific effects
matrix[N*(K-1),T] Z; # data for alt constant effects
}
transformed data {
}
parameters {
matrix[D, K-1] beta; # individual specific coefs
matrix[G, K] gamma; # choice specific coefs for alt-specific variables
vector[T] theta; # choice constant coefs for alt-specific variables
}
transformed parameters{
}
model {
matrix[N, K-1] Vx; # Utility for individual vars
vector[N*K] Vy0;
matrix[N, K-1] Vy; # Utility for alt-specific/alt-varying vars
vector[N*(K-1)] Vz0;
matrix[N, (K-1)] Vz; # Utility for alt-specific/alt-constant vars
matrix[N,K-1] V; # combined utilities
vector[N] baseProbVec; # reference group probabilities
real ll0; # intermediate log likelihood
real loglik; # final log likelihood
to_vector(beta) ~ normal(0, 10); # diffuse priors on coefficients
to_vector(gamma) ~ normal(0, 10); # diffuse priors on coefficients
to_vector(theta) ~ normal(0, 10); # diffuse priors on coefficients
Vx = X * beta;
for(alt in 1:K){
vector[G] par;
int start;
int end;
par = gamma[,alt];
start = N*alt-N+1;
end = N*alt;
Vy0[start:end] = Y[start:end,] * par;
if(alt>1) Vy[,alt-1] = Vy0[start:end] - Vy0[1:N];
}
Vz0 = Z * theta;
for(alt in 1:(K-1)){
int start;
int end;
start = N*alt-N+1;
end = N*alt;
Vz[,alt] = Vz0[start:end];
}
V = Vx + Vy + Vz;
for(n in 1:N) baseProbVec[n] = 1/(1 + sum(exp(V[n])));
ll0 = dot_product(to_vector(V), choice[(N+1):(N*K)]); # just going to assume no neg index
loglik = sum(log(baseProbVec)) + ll0;
target += loglik;
}
generated quantities {
matrix[N,K-1] fitted_nonref;
vector[N] fitted_ref;
matrix[N,K] fitted;
matrix[N,K-1] Vx; # Utility for individual vars
vector[N*K] Vy0;
matrix[N,K-1] Vy; # Utility for alt-specific/alt-varying vars
vector[N*(K-1)] Vz0;
matrix[N, (K-1)] Vz; # Utility for alt-specific/alt-constant vars
matrix[N,K-1] V; # combined utilities
vector[N] baseProbVec; # reference group probabilities
Vx = X * beta;
for(alt in 1:K){
vector[G] par;
int start;
int end;
par = gamma[,alt];
start = N*alt-N+1;
end = N*alt;
Vy0[start:end] = Y[start:end,] * par;
if(alt>1) Vy[,alt-1] = Vy0[start:end] - Vy0[1:N];
}
Vz0 = Z * theta;
for(alt in 1:(K-1)){
int start;
int end;
start = N*alt-N+1;
end = N*alt;
Vz[,alt] = Vz0[start:end];
}
V = Vx + Vy + Vz;
for(n in 1:N) baseProbVec[n] = 1/(1 + sum(exp(V[n])));
fitted_nonref = exp(V) .* rep_matrix(baseProbVec, K-1);
for(n in 1:N) fitted_ref[n] = 1-sum(fitted_nonref[n]);
fitted = append_col(fitted_ref, fitted_nonref);
}