teratail
回答率
85.47%
質問するログイン新規登録
トップに関する質問rstanを利用してベイズ推定を行いたい。

編集履歴

質問編集履歴

1

自分で書いたコード、参考にしたコードとその参考サイトのURLを記載しました。

2019/01/22 06:48

投稿

kzzzzz
kzzzzz

スコア10

test CHANGED
File without changes
test CHANGED
@@ -37,3 +37,247 @@
37
37
 
38
38
 
39
39
  ということを知りたいです。
40
+
41
+
42
+
43
+ ### 参考にしたコード
44
+
45
+
46
+
47
+ model.stan
48
+
49
+ ```
50
+
51
+ data {
52
+
53
+ int t;
54
+
55
+ int Dat[t];
56
+
57
+ }
58
+
59
+
60
+
61
+ parameters {
62
+
63
+ vector <lower = -5, upper = 5>[t] mu;
64
+
65
+ real<lower=0> sigma;
66
+
67
+ }
68
+
69
+
70
+
71
+ transformed parameters {
72
+
73
+ vector<lower=0, upper=1>[t] theta;
74
+
75
+ theta = inv_logit(mu);
76
+
77
+ }
78
+
79
+
80
+
81
+ model {
82
+
83
+ sigma ~ normal(0, 1);
84
+
85
+ for (i in 2:t){
86
+
87
+ mu[i] ~ normal(mu[i-1], sigma);
88
+
89
+ }
90
+
91
+ Dat ~ bernoulli(theta);
92
+
93
+ }
94
+
95
+
96
+
97
+ ```
98
+
99
+
100
+
101
+ Rコード
102
+
103
+ ```
104
+
105
+ #data to be passed on stan
106
+
107
+ t=length(Dat)
108
+
109
+ datastan=list(t=t,Dat=Dat)
110
+
111
+
112
+
113
+ #
114
+
115
+ # stan model
116
+
117
+ #
118
+
119
+ #model<-stan_model("GaussianLogitDynamic.stan")
120
+
121
+ #saveRDS(model,"GaussianLogitDynamic.rds")
122
+
123
+ GLD<-readRDS("GaussianLogitDynamic.rds")
124
+
125
+
126
+
127
+ #
128
+
129
+ # stan fit
130
+
131
+ #
132
+
133
+ fit = sampling(GLD, data = datastan, chains = 4, iter = 4000, warmup = 2000, thin=4, seed=1234)
134
+
135
+
136
+
137
+ #Gelman-Rubin
138
+
139
+ launch_shinystan(fit)
140
+
141
+ ```
142
+
143
+ URL
144
+
145
+ ```
146
+
147
+ https://mrunadon.github.io/GaussianLogitDynamic/
148
+
149
+ ```
150
+
151
+
152
+
153
+ ### 自分の書いたコード
154
+
155
+ pre_model.stan
156
+
157
+ ```
158
+
159
+ data{
160
+
161
+ int t;
162
+
163
+ int Dat[t];
164
+
165
+ }
166
+
167
+
168
+
169
+ parameters{
170
+
171
+ real v1;
172
+
173
+ real v2;
174
+
175
+ real<lower=0> sigma;
176
+
177
+ }
178
+
179
+
180
+
181
+ transformed parameters{
182
+
183
+ vector<lower=0,upper=1>[t] x;
184
+
185
+ vector<lower=0,upper=1> [t] y;
186
+
187
+ vector<lower=0,upper=1>[t] theta;
188
+
189
+ x[1]=0;
190
+
191
+ for(i in 2:t){
192
+
193
+ if(Dat[i]==1){
194
+
195
+ x[i]=x[i-1]+v2;
196
+
197
+ }
198
+
199
+ if(Dat[i]==0){
200
+
201
+ x[i]=x[i-1]+v1;
202
+
203
+ }
204
+
205
+ }
206
+
207
+
208
+
209
+ y[1]=0.5;
210
+
211
+ theta[t]=(1-y[t])*x[t];
212
+
213
+ }
214
+
215
+
216
+
217
+ model{
218
+
219
+ v1 ~ normal(0,1);
220
+
221
+ v2 ~ normal(0,1);
222
+
223
+ sigma ~ normal(0,1);
224
+
225
+ for(i in 2:t){
226
+
227
+ y[i] ~ normal(y[i-1],sigma);
228
+
229
+ }
230
+
231
+ Dat ~ bernoulli(theta);
232
+
233
+ }
234
+
235
+
236
+
237
+ ```
238
+
239
+
240
+
241
+ Rコード
242
+
243
+ データは事前に指定しています。
244
+
245
+ ```
246
+
247
+ #data to be passed on stan
248
+
249
+ t=length(Dat)
250
+
251
+ datastan=list(t=t,Dat=Dat)
252
+
253
+
254
+
255
+ #
256
+
257
+ # stan model
258
+
259
+ #
260
+
261
+ #model<-stan_model("pre_model.stan")
262
+
263
+ #saveRDS(model,"pre_model.rds")
264
+
265
+ GLD<-readRDS("pre_model.rds")
266
+
267
+
268
+
269
+ #
270
+
271
+ # stan fit
272
+
273
+ #
274
+
275
+ fit = sampling(GLD, data = datastan, chains = 4, iter = 4000, warmup = 2000, thin=4, seed=1234)
276
+
277
+
278
+
279
+ #Gelman-Rubin
280
+
281
+ launch_shinystan(fit)
282
+
283
+ ```