1 ////////////////////////////////////////////////////////////////////
2 // cMCMCprobitChange.cc is C++ code to estimate a probit 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 // Written 07/06/2007
11 // Modifieed 11/02/2009
12 // Included 1/21/2011
13 //////////////////////////////////////////////////////////////////////////
14
15 #ifndef CMCMCPROBITCHANGE_CC
16 #define CMCMCPROBITCHANGE_CC
17
18 #include "MCMCrng.h"
19 #include "MCMCfcds.h"
20 #include "matrix.h"
21 #include "distributions.h"
22 #include "stat.h"
23 #include "la.h"
24 #include "ide.h"
25 #include "smath.h"
26 #include "rng.h"
27 #include "mersenne.h"
28 #include "lecuyer.h"
29
30 #include <R.h>
31 #include <R_ext/Utils.h>
32
33
34 using namespace std;
35 using namespace scythe;
36
37 // probit state sampler
38 template <typename RNGTYPE>
39
probit_state_sampler(rng<RNGTYPE> & stream,const int m,const Matrix<double> & Y,const Matrix<double> & X,const Matrix<double> & beta,const Matrix<double> & P)40 Matrix<> probit_state_sampler(rng<RNGTYPE>& stream,
41 const int m,
42 const Matrix<double>& Y,
43 const Matrix<double>& X,
44 const Matrix<double>& beta,
45 const Matrix<double>& P){
46
47 const int ns = m + 1;
48 const int n = Y.rows();
49 Matrix<double> F = Matrix<double>(n, ns);
50 Matrix<double> pr1 = Matrix<double>(ns, 1);
51 pr1[0] = 1;
52 Matrix<double> py(ns, 1);
53 Matrix<double> pstyt1(ns, 1);
54
55 for (int t=0; t<n ; ++t){
56 int yt = (int) Y[t];
57 Matrix<double> mu = X(t,_)*::t(beta);
58 for (int j=0; j<ns ; ++j){
59 py[j] = dbinom(yt, 1, pnorm(mu[j], 0, 1));
60 }
61 if (t==0) pstyt1 = pr1;
62 else {
63 pstyt1 = ::t(F(t-1,_)*P);
64 }
65 Matrix<double> unnorm_pstyt = pstyt1%py;
66 const Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
67 for (int j=0; j<ns ; ++j){
68 F(t,j) = pstyt(j);
69 }
70 }
71
72 Matrix<int> s(n, 1);
73 Matrix<double> ps = Matrix<double>(n, ns);
74 ps(n-1,_) = F(n-1,_);
75 s(n-1) = ns;
76
77 Matrix<double> pstyn = Matrix<double>(ns, 1);
78 double pone = 0.0;
79 int t = n-2;
80 while (t >= 0){
81 int st = s(t+1);
82 Matrix<double> Pst_1 = ::t(P(_,st-1));
83 Matrix<double> unnorm_pstyn = F(t,_)%Pst_1;
84 pstyn = unnorm_pstyn/sum(unnorm_pstyn);
85 if (st==1) s(t) = 1;
86 else{
87 pone = pstyn(st-2);
88 if(stream.runif() < pone) s(t) = st-1;
89 else s(t) = st;
90 }
91 ps(t,_) = pstyn;
92 --t;
93 }
94
95 Matrix<double> Sout(n, ns+1);
96 Sout(_, 0) = s(_,0);
97 for (int j = 0; j<ns; ++j){
98 Sout(_,j+1) = ps(_, j);
99 }
100
101 return Sout;
102
103 }
104
105 ////////////////////////////////////////////
106 // cMCMCprobitChangepoint implementation.
107 ////////////////////////////////////////////
108 template <typename RNGTYPE>
MCMCprobitChange_impl(rng<RNGTYPE> & stream,const int m,const Matrix<> & Y,const Matrix<> & X,Matrix<> & beta,Matrix<> & P,Matrix<> & b0,Matrix<> & B0,const Matrix<> & A0,int burnin,int mcmc,int thin,int verbose,bool chib,Matrix<> & beta_store,Matrix<> & Z_store,Matrix<> & P_store,Matrix<> & ps_store,Matrix<int> & s_store,double & logmarglike,double & loglike)109 void MCMCprobitChange_impl(rng<RNGTYPE>& stream,
110 const int m,
111 const Matrix<>& Y, const Matrix<>& X,
112 Matrix<>& beta, Matrix<>& P,
113 Matrix<>& b0, Matrix<>& B0,
114 const Matrix<>& A0,
115 int burnin, int mcmc, int thin,
116 int verbose, bool chib,
117 Matrix<>& beta_store, Matrix<>& Z_store,
118 Matrix<>& P_store, Matrix<>& ps_store, Matrix<int>& s_store,
119 double& logmarglike, double& loglike)
120 {
121 const int tot_iter = burnin + mcmc;
122 const int nstore = mcmc / thin;
123 const int n = Y.rows();
124 const int ns = m + 1;
125 const int k = X.cols();
126 const Matrix<> B0inv = invpd(B0);
127 Matrix<> Z(n, 1);
128
129 int count = 0;
130 for (int iter = 0; iter < tot_iter; ++iter){
131
132 // 1. Sample s
133 Matrix<> Sout = probit_state_sampler(stream, m, Y, X, beta, P);
134 Matrix<int> s = Sout(_, 0);
135 Matrix <double> ps(n, ns);
136 for (int j = 0; j<ns; ++j){
137 ps(_,j) = Sout(_,j+1);
138 }
139
140 // 2. Sample Z
141 for ( int i = 0; i < n; ++i) {
142 const Matrix<> mu = X(i,_)*t(beta(s[i]-1,_));
143 double muj = mu[0];
144 if(muj>200){
145 muj = 200;
146 }
147 if(muj<-200){
148 muj = -200;
149 }
150 if (Y[i] == 1.0)
151 Z[i] = stream.rtbnorm_combo(muj, 1.0, 0);
152 if (Y[i] == 0.0)
153 Z[i] = stream.rtanorm_combo(muj, 1.0, 0);
154 }
155
156 // 3. Sample beta
157 int beta_count = 0;
158 Matrix<int> nstate(ns, 1);
159 for (int j = 0; j <ns ; ++j){
160 for (int i = 0; i<n; ++i){
161 if (s[i] == (j+1)) {
162 nstate[j] = nstate[j] + 1;
163 }
164 }
165 beta_count = beta_count + nstate[j];
166 const Matrix<double> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
167 const Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
168 const Matrix<double> XpX = t(Xj)*Xj;
169 const Matrix<double> XpZ = t(Xj)*Zj;
170 beta(j,_) = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
171 }
172
173 // 4. Sample P
174 double shape1 = 0;
175 double shape2 = 0;
176 P(ns-1, ns-1) = 1;
177
178 for (int j =0; j< (ns-1); ++j){
179 shape1 = A0(j,j) + nstate[j] - 1;
180 shape2 = A0(j,j+1) + 1; //
181 P(j,j) = stream.rbeta(shape1, shape2);
182 P(j,j+1) = 1 - P(j,j);
183 }
184
185 // load draws into sample array
186 if (iter >= burnin && ((iter % thin)==0)){
187 Matrix<double> tbeta = ::t(beta);
188 for (int i=0; i<(ns*k); ++i)
189 beta_store(count,i) = tbeta[i];
190 for (int j=0; j<ns*ns; ++j)
191 P_store(count,j)= P[j];
192 for (int l=0; l<n ; ++l)
193 ps_store(l,_) = ps_store(l,_) + ps(l,_);
194 s_store(count,_) = s;
195 Z_store(count,_) = Z;
196
197 ++count;
198
199 } // end of if(iter...)
200
201
202 // print output to stdout
203 if(verbose > 0 && iter % verbose == 0){
204 Rprintf("\n\nMCMCprobitChange iteration %i of %i \n", (iter+1), tot_iter);
205 for (int j = 0;j<ns; ++j){
206 Rprintf("\n The number of observations in state %i is %10.5d", j+1, nstate[j]);
207 }
208 for (int j = 0; j<ns; ++j){
209 Rprintf("\n beta %i is ", j);
210 for (int i = 0; i<k; ++i){
211 Rprintf("%10.5f", beta(j, i));
212 }
213 }
214 }
215
216 R_CheckUserInterrupt();
217
218 }// end MCMC loop
219
220 if(chib==1){
221 Matrix<double> betast = meanc(beta_store);
222 Matrix<double, Row> beta_st(ns, k);
223 for (int j = 0; j<ns*k; ++j){
224 beta_st[j] = betast[j];
225 }
226 Matrix<double> P_vec_st = meanc(P_store);
227 const Matrix<double> P_st(ns, ns);
228 for (int j = 0; j< ns*ns; ++j){
229 P_st[j] = P_vec_st[j];
230 }
231
232 // 1. beta
233 Matrix<double> density_beta(nstore, ns);
234 for (int iter = 0; iter<nstore; ++iter){
235 Matrix<int> nstate(ns, 1);
236 int beta_count = 0;
237 const Matrix<double> Z(n, 1);
238 Z(_,0) = Z_store(iter,_);
239 for (int j = 0; j<ns ; ++j){
240 for (int i = 0; i<n; ++i){
241 if (s_store(iter, i) == (j+1)) {
242 nstate[j] = nstate[j] + 1;
243 }
244 }
245 beta_count = beta_count + nstate[j];
246 const Matrix<double> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
247 const Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
248 const Matrix<double> XpX = (::t(Xj)*Xj);
249 const Matrix<double> XpZ = (::t(Xj)*Zj);
250 const Matrix<double> Bn = invpd(B0 + XpX);
251 const Matrix<double> bn = Bn*gaxpy(B0, b0, XpZ);
252 density_beta(iter, j) = exp(lndmvn(::t(beta_st(j,_)), bn, Bn));
253 }
254 }
255
256 double pdf_beta = log(prod(meanc(density_beta)));
257
258 // 2. P
259 Matrix<double> density_P(nstore, ns);
260 for (int iter = 0; iter < nstore; ++iter){
261 Matrix <double> Sout = probit_state_sampler(stream, m, Y, X, beta_st, P);
262 Matrix <double> s = Sout(_, 0);
263 Matrix <double> ps(n, ns);
264 for (int j = 0; j<ns; ++j){
265 ps(_,j) = Sout(_,j+1);
266 }
267
268 double shape1 = 0;
269 double shape2 = 0;
270 P(ns-1, ns-1) = 1;
271 Matrix <double> P_addN(ns, 1);
272 for (int j = 0; j<ns; ++j){
273 for (int i = 0; i<n; ++i){
274 if (s[i] == (j+1)) {
275 P_addN[j] = P_addN[j] + 1;
276 }
277 }
278 }
279 for (int j =0; j< (ns-1); ++j){
280 shape1 = A0(j,j) + P_addN[j] - 1;
281 shape2 = A0(j,j+1) + 1; //
282 P(j,j) = stream.rbeta(shape1, shape2);
283 P(j,j+1) = 1 - P(j,j);
284 density_P(iter, j) = dbeta(P_st(j,j), shape1, shape2);
285 }
286 density_P(iter, ns-1) = 1;
287
288 }
289 double pdf_P = log(prod(meanc(density_P)));
290
291 // likelihood
292 Matrix<double> F = Matrix<double>(n, ns);
293 Matrix<double> like(n, 1);
294 Matrix<double> pr1 = Matrix<double>(ns, 1);
295 pr1[0] = 1;
296 Matrix<double> py(ns, 1);
297 Matrix<double> pstyt1(ns, 1);
298
299 for (int t=0; t<n ; ++t){
300 int yt = (int) Y[t];
301 Matrix<double> mu = X(t,_)*::t(beta_st);
302 for (int j=0; j<ns ; ++j){
303 py[j] = dbinom(yt, 1, pnorm(mu[j], 0, 1));
304 }
305 if (t==0) pstyt1 = pr1;
306 else {
307 pstyt1 = ::t(F(t-1,_)*P_st);
308 }
309 Matrix<double> unnorm_pstyt = pstyt1%py;
310 Matrix<double> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
311 for (int j=0; j<ns ; ++j){
312 F(t,j) = pstyt(j);
313 }
314 like[t] = sum(unnorm_pstyt);
315 }
316
317 loglike = sum(log(like));
318
319 //////////////////////
320 // log prior ordinate
321 //////////////////////
322 Matrix<double> density_beta_prior(ns, 1);
323 Matrix<double> density_P_prior(ns, 1);
324 density_P[ns-1] = 1; //
325
326 for (int j=0; j<ns ; ++j){
327 density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
328 }
329
330 for (int j =0; j< (ns-1); ++j){
331 density_P_prior[j] = scythe::lndbeta1(P_st(j,j), A0(j,j), A0(j,j+1));
332 }
333
334 // compute marginal likelihood
335 double logprior = sum(density_beta_prior) + sum(density_P_prior);
336 logmarglike = (loglike + logprior) - (pdf_beta + pdf_P);
337
338 if (verbose >0 ){
339 Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
340 Rprintf("loglike = %10.5f\n", loglike);
341 Rprintf("log_prior = %10.5f\n", logprior);
342 Rprintf("log_beta = %10.5f\n", pdf_beta);
343 Rprintf("log_P = %10.5f\n", pdf_P);
344 }
345 } // end of marginal likelihood
346 }//end extern "C"
347
348 extern "C"{
cMCMCprobitChange(double * betaout,double * Pout,double * psout,double * sout,const double * Ydata,const int * Yrow,const int * Ycol,const double * Xdata,const int * Xrow,const int * Xcol,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 * betastart,const double * Pstart,const double * a,const double * b,const double * b0data,const double * B0data,const double * A0data,double * logmarglikeholder,double * loglikeholder,const int * chib)349 void cMCMCprobitChange(double *betaout, double *Pout, double *psout, double *sout,
350 const double *Ydata, const int *Yrow, const int *Ycol,
351 const double *Xdata, const int *Xrow, const int *Xcol,
352 const int *m,
353 const int *burnin, const int *mcmc, const int *thin, const int *verbose,
354 const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
355 const double *betastart, const double *Pstart,
356 const double *a, const double *b,
357 const double *b0data, const double *B0data,
358 const double *A0data,
359 double *logmarglikeholder, double *loglikeholder,
360 const int *chib){
361
362 // pull together Matrix objects
363 const Matrix <> Y(*Yrow, *Ycol, Ydata);
364 const Matrix <> X(*Xrow, *Xcol, Xdata);
365 const int nstore = *mcmc / *thin;
366 const int n = Y.rows();
367 const int k = X.cols();
368 const int ns = *m + 1;
369
370 // generate starting values
371 Matrix <> beta(ns, k, betastart);
372 Matrix <> P(ns, ns, Pstart);
373 Matrix <> b0(k, 1, b0data);
374 Matrix <> B0(k, k, B0data);
375 const Matrix <> A0(ns, ns, A0data);
376 double logmarglike;
377 double loglike;
378
379 // storage matrices
380 Matrix<> beta_store(nstore, ns*k);
381 Matrix<> Z_store(nstore, n);
382 Matrix<> P_store(nstore, ns*ns);
383 Matrix<> ps_store(n, ns);
384 Matrix<int> s_store(nstore, n);
385
386 MCMCPACK_PASSRNG2MODEL(MCMCprobitChange_impl, *m,
387 Y, X, beta, P, b0, B0, A0,
388 *burnin, *mcmc, *thin, *verbose, *chib,
389 beta_store, Z_store,
390 P_store, ps_store, s_store,
391 logmarglike, loglike);
392 logmarglikeholder[0] = logmarglike;
393 loglikeholder[0] = loglike;
394
395 // return output
396 for (int i = 0; i<(nstore*ns*k); ++i){
397 betaout[i] = beta_store[i];
398 }
399 for (int i = 0; i<(nstore*ns*ns); ++i){
400 Pout[i] = P_store[i];
401 }
402 for (int i = 0; i<(n*ns); ++i){
403 psout[i] = ps_store[i];
404 }
405 for (int i = 0; i<(nstore*n); ++i){
406 sout[i] = s_store[i];
407 }
408
409 }
410 }
411 #endif
412
413
414