1 ////////////////////////////////////////////////////////////////////
2 // cMCMCbinaryChange.cc is C++ code to estimate a binary changepoint model
3 // with a beta prior
4 //
5 // Jong Hee Park
6 // Department of Political Science and International Relations
7 // Seoul National University
8 // jongheepark@snu.ac.kr
9 //
10 // 03/03/2009 Written
11 ////////////////////////////////////////////////////////////////////
12
13 #ifndef CMCMCBINARYCHANGE_CC
14 #define CMCMCBINARYCHANGE_CC
15
16
17 #include "MCMCrng.h"
18 #include "MCMCfcds.h"
19 #include "matrix.h"
20 #include "distributions.h"
21 #include "stat.h"
22 #include "la.h"
23 #include "ide.h"
24 #include "smath.h"
25 #include "rng.h"
26 #include "mersenne.h"
27 #include "lecuyer.h"
28
29 #include <R.h>
30 #include <R_ext/Utils.h>
31
32
33 template <typename RNGTYPE>
34 // bianry state sampler
binary_state_sampler(rng<RNGTYPE> & stream,const int m,const Matrix<> & Y,const Matrix<> & theta,const Matrix<> & P)35 Matrix<> binary_state_sampler(rng<RNGTYPE>& stream,
36 const int m,
37 const Matrix<>& Y,
38 const Matrix<>& theta,
39 const Matrix<>& P){
40
41 const int ns = m + 1;
42 const int n = Y.rows();
43 Matrix<> F = Matrix<>(n, ns);
44 Matrix<> pr1 = Matrix<>(ns, 1);
45 pr1[0] = 1;
46 Matrix<> py(ns, 1);
47 Matrix<> pstyt1(ns, 1);
48
49 for (int t=0; t<n ; ++t){
50 int yt = (int) Y[t];
51 for (int j=0; j<ns ; ++j){
52 py[j] = dbinom(yt, 1, theta[j]);
53 }
54 if (t==0) pstyt1 = pr1;
55 else {
56 pstyt1 = ::t(F(t-1,_)*P);
57 }
58 Matrix<> unnorm_pstyt = pstyt1%py;
59 const Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
60 for (int j=0; j<ns ; ++j){
61 F(t,j) = pstyt(j);
62 }
63 }
64
65 Matrix<int> s(n, 1);
66 Matrix<> ps = Matrix<>(n, ns);
67 ps(n-1,_) = F(n-1,_);
68 s(n-1) = ns;
69
70 Matrix<> pstyn = Matrix<>(ns, 1);
71 double pone = 0.0;
72 int t = n-2;
73 while (t >= 0){
74 int st = s(t+1);
75 Matrix<> Pst_1 = ::t(P(_,st-1));
76 Matrix<> unnorm_pstyn = F(t,_)%Pst_1;
77 pstyn = unnorm_pstyn/sum(unnorm_pstyn);
78 if (st==1) s(t) = 1;
79 else{
80 pone = pstyn(st-2);
81 if(stream.runif() < pone) s(t) = st-1;
82 else s(t) = st;
83 }
84 ps(t,_) = pstyn;
85 --t;
86 }
87
88 Matrix<> Sout(n, ns+1);
89 Sout(_, 0) = s(_,0);
90 for (int j = 0; j<ns; ++j){
91 Sout(_,j+1) = ps(_, j);
92 }
93
94 return Sout;
95
96 }
97
98 ////////////////////////////////////////////
99 // cMCMCbinaryChangepoint implementation.
100 ////////////////////////////////////////////
101 template <typename RNGTYPE>
MCMCbinaryChange_impl(rng<RNGTYPE> & stream,const Matrix<> & Y,Matrix<> & phi,Matrix<> & P,const Matrix<> & A0,const double m,const double c0,const double d0,int burnin,int mcmc,int thin,int verbose,bool chib,Matrix<> & phi_store,Matrix<> & P_store,Matrix<> & ps_store,Matrix<> & s_store,double & logmarglike)102 void MCMCbinaryChange_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
103 Matrix<>& phi, Matrix<>& P, const Matrix<>& A0,
104 const double m, const double c0, const double d0,
105 int burnin, int mcmc, int thin,
106 int verbose, bool chib,
107 Matrix<>& phi_store, Matrix<>& P_store,
108 Matrix<>& ps_store, Matrix<>& s_store,
109 double& logmarglike)
110 {
111 const int tot_iter = burnin + mcmc;
112 const int nstore = mcmc / thin;
113 const int n = Y.rows();
114 const int ns = m + 1;
115
116 //MCMC loop
117 int count = 0;
118 for (int iter = 0; iter < tot_iter; ++iter){
119
120 //////////////////////
121 // 1. Sample s
122 //////////////////////
123 Matrix<> Sout = binary_state_sampler(stream, m, Y, phi, P);
124 Matrix<> s = Sout(_, 0);
125 Matrix<> ps(n, ns);
126 for (int j = 0; j<ns; ++j){
127 ps(_,j) = Sout(_,j+1);
128 }
129
130 //////////////////////
131 // 2. Sample phi
132 //////////////////////
133 Matrix<> addY(ns, 1);
134 Matrix<> addN(ns, 1);
135
136 for (int j = 0; j<ns; ++j){
137 for (int i = 0; i<n; ++i){
138 if (s[i] == (j+1)) {
139 addN[j] = addN[j] + 1;
140 addY[j] = addY[j] + Y[i];
141 }
142 }
143 double c1 = addY[j] + c0;
144 double d1 = addN[j] - addY[j] + d0;
145 phi[j] = stream.rbeta(c1, d1);
146 }
147
148 //////////////////////
149 // 3. Sample P
150 //////////////////////
151 double shape1 = 0;
152 double shape2 = 0;
153 P(ns-1, ns-1) = 1;
154
155 for (int j =0; j< (ns-1); ++j){
156 shape1 = A0(j,j) + addN[j] - 1;
157 shape2 = A0(j,j+1) + 1;
158 P(j,j) = stream.rbeta(shape1, shape2);
159 P(j,j+1) = 1 - P(j,j);
160 }
161
162 if (iter >= burnin && ((iter % thin)==0)){
163 for (int i=0; i<ns; ++i)
164 phi_store(count,i) = phi[i];
165 for (int j=0; j<ns*ns; ++j)
166 P_store(count,j)= P[j];
167 s_store(count,_) = s;
168 for (int l=0; l<n ; ++l)
169 ps_store(l,_) = ps_store(l,_) + ps(l,_);
170
171 ++count;
172
173 }
174 if(verbose > 0 && iter % verbose == 0){
175 Rprintf("\n\n MCMCbinaryChange iteration %i of %i", (iter+1), tot_iter);
176 for (int j = 0;j<ns; ++j){
177 Rprintf("\n The number of observations in state %i is %10.5f", j + 1, addN[j]);
178 }
179 for (int j=0; j<ns; ++j)
180 Rprintf("\n phi in state %i is %10.5f", j + 1, phi[j]);
181 }
182
183 R_CheckUserInterrupt();
184
185 }// end MCMC loop
186
187 if(chib ==1){
188
189 Matrix<> phi_st = meanc(phi_store);
190 Matrix<> P_vec_st = meanc(P_store);
191 const Matrix<> P_st(ns, ns);
192 for (int j = 0; j< ns*ns; ++j){
193 P_st[j] = P_vec_st[j];
194 }
195
196 //////////////////////
197 // phi
198 //////////////////////
199 Matrix<> density_phi(nstore, ns);
200 // Matrix<> density_phi_j(ns, 1);
201 for (int iter = 0; iter<nstore; ++iter){
202
203 Matrix<> addY(ns, 1);
204 Matrix<> addN(ns, 1);
205
206 for (int j = 0; j<ns; ++j){
207 for (int i = 0; i<n; ++i){
208 if (s_store(iter, i) == (j+1)) {
209 addN[j] = addN[j] + 1;
210 addY[j] = addY[j] + Y[i];
211 }
212 }
213 double c1 = addY[j] + c0;
214 double d1 = addN[j] - addY[j] + d0;
215 density_phi(iter, j) = dbeta(phi_st[j], c1, d1);
216 }
217 }
218 double pdf_phi = log(prod(meanc(density_phi)));
219
220 //////////////////////
221 // P
222 //////////////////////
223 Matrix<> density_P(nstore, ns);
224 for (int iter = 0; iter < nstore; ++iter){
225 Matrix<> Sout = binary_state_sampler(stream, m, Y, phi_st, P);
226 Matrix <> s = Sout(_, 0);
227 Matrix <> ps(n, ns);
228 for (int j = 0; j<ns; ++j){
229 ps(_,j) = Sout(_,j+1);
230 }
231
232 double shape1 = 0;
233 double shape2 = 0;
234 P(ns-1, ns-1) = 1;
235 Matrix <> P_addN(ns, 1);
236 for (int j = 0; j<ns; ++j){
237 for (int i = 0; i<n; ++i){
238 if (s[i] == (j+1)) {
239 P_addN[j] = P_addN[j] + 1;
240 }
241 }
242 }
243
244 for (int j =0; j< (ns-1); ++j){
245 shape1 = A0(j,j) + P_addN[j] - 1;
246 shape2 = A0(j,j+1) + 1;
247 P(j,j) = stream.rbeta(shape1, shape2);
248 P(j,j+1) = 1 - P(j,j);
249 density_P(iter, j) = dbeta(P_st(j,j), shape1, shape2);
250 }
251 density_P(iter, ns-1) = 1;
252
253 }
254 double pdf_P = log(prod(meanc(density_P)));
255
256 //////////////////////
257 // likelihood
258 //////////////////////
259 Matrix<> F(n, ns);
260 Matrix<> like(n, 1);
261 Matrix<> pr1(ns, 1);
262 pr1[0] = 1;
263 Matrix<> py(ns, 1);
264 Matrix<> pstyt1(ns, 1);
265
266 for (int t=0; t<n ; ++t){
267 int yt = (int) Y[t];
268 for (int j=0; j<ns ; ++j){
269 py[j] = dbinom(yt, 1, phi_st[j]);
270 }
271 if (t==0) pstyt1 = pr1;
272 else {
273 pstyt1 = ::t(F(t-1,_)*P_st);
274 }
275 Matrix<> unnorm_pstyt = pstyt1%py;
276 Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
277 for (int j=0; j<ns ; ++j){
278 F(t,j) = pstyt(j);
279 }
280 like[t] = sum(unnorm_pstyt);
281 }
282
283 double loglike = sum(log(like));
284
285 //////////////////////
286 // prior
287 //////////////////////
288 Matrix<> density_phi_prior(ns, 1);
289 Matrix<> density_P_prior(ns, 1);
290
291 for (int j=0; j<ns ; ++j){
292 density_phi_prior[j] = scythe::lndbeta1(phi_st[j], c0, d0);
293 }
294
295 for (int j =0; j< (ns-1); ++j){
296 density_P_prior[j] = scythe::lndbeta1(P_st(j,j), A0(j,j), A0(j,j+1));
297 }
298
299 // compute marginal likelihood
300 double logprior = sum(density_phi_prior) + sum(density_P_prior);
301 logmarglike = (loglike + logprior) - (pdf_phi + pdf_P);
302 if(verbose > 0){
303 Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
304 Rprintf("loglike = %10.5f\n", loglike);
305 Rprintf("log_prior = %10.5f\n", logprior);
306 Rprintf("log_phi = %10.5f\n", pdf_phi);
307 Rprintf("log_P = %10.5f\n", pdf_P);
308 }
309
310 }
311 }
312 ////////////////////////////////////////////
313 // Start cMCMCbinaryChangepoint function
314 ///////////////////////////////////////////
315 extern "C"{
cMCMCbinaryChange(double * phiout,double * Pout,double * psout,double * sout,const double * Ydata,const int * Yrow,const int * Ycol,const int * m,const int * burnin,const int * mcmc,const int * thin,const int * verbose,const int * uselecuyer,const int * seedarray,const int * lecuyerstream,const double * phistart,const double * Pstart,const double * a,const double * b,const double * c0,const double * d0,const double * A0data,double * logmarglikeholder,const int * chib)316 void cMCMCbinaryChange(double *phiout, double *Pout, double *psout, double *sout,
317 const double *Ydata, const int *Yrow, const int *Ycol, const int *m,
318 const int *burnin, const int *mcmc, const int *thin, const int *verbose,
319 const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
320 const double *phistart, const double *Pstart,
321 const double *a, const double *b, const double *c0, const double *d0,
322 const double *A0data, double *logmarglikeholder,
323 const int *chib){
324
325 // pull together Matrix objects
326 const Matrix <> Y(*Yrow, *Ycol, Ydata);
327 const int nstore = *mcmc / *thin;
328 const int n = Y.rows();
329 const int ns = *m + 1;
330
331 // generate starting values
332 Matrix <> phi(ns, 1, phistart);
333 const Matrix <> A0(ns, ns, A0data);
334 Matrix <> P(ns, ns, Pstart);
335
336 double logmarglike;
337
338 // storage matrices
339 Matrix<> phi_store(nstore, ns);
340 Matrix<> P_store(nstore, ns*ns);
341 Matrix<> ps_store(n, ns);
342 Matrix<> s_store(nstore, n);
343
344 MCMCPACK_PASSRNG2MODEL(MCMCbinaryChange_impl, Y,
345 phi, P, A0, *m, *c0, *d0,
346 *burnin, *mcmc, *thin, *verbose, *chib,
347 phi_store, P_store, ps_store, s_store,
348 logmarglike)
349
350 logmarglikeholder[0] = logmarglike;
351
352 // return output
353 for (int i = 0; i<(nstore*ns); ++i){
354 phiout[i] = phi_store[i];
355 }
356 for (int i = 0; i<(nstore*ns*ns); ++i){
357 Pout[i] = P_store[i];
358 }
359 for (int i = 0; i<(n*ns); ++i){
360 psout[i] = ps_store[i];
361 }
362 for (int i = 0; i<(nstore*n); ++i){
363 sout[i] = s_store[i];
364 }
365
366 }
367 }
368
369 #endif
370
371