1 ////////////////////////////////////////////////////////////////////
2 // cMCMCoprobitChange.cc is C++ code to estimate a oprobit changepoint model
3 // with linear approximation
4 //
5 // Jong Hee Park
6 // Department of Political Science and International Relations
7 // Seoul National University
8 // jongheepark@snu.ac.kr
9 //
10 // 07/06/2007 Written
11 // 11/02/2009 Modified
12 ////////////////////////////////////////////////////////////////////
13
14 #ifndef CMCMCOPROBITCHANGE_CC
15 #define CMCMCOPROBITCHANGE_CC
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 using namespace std;
34 using namespace scythe;
35
36 // density function for truncated normal
dtnormLX(const double x,const double mean,const double sd,const double lower,const double upper)37 static double dtnormLX(const double x,
38 const double mean,
39 const double sd,
40 const double lower,
41 const double upper){
42 double out = 0.0;
43 if (x>lower && x<upper){
44 double numer = dnorm(x, mean, sd);
45 double denom = pnorm(upper, mean, sd) - pnorm(lower, mean, sd);
46 out = numer/denom;
47 }
48 else{
49 Rprintf("\n x input for dtnormLX() %10.5f is out of bounds %10.5f %10.5f ", x, lower, upper, "\n");
50 out = 1;
51 }
52 return out;
53 }
54
55 // likelihood of oprobit
oprobit_pdfLX(const int ncat,const int Y,double Xbeta,const Matrix<> & gamma)56 static double oprobit_pdfLX(const int ncat,
57 const int Y,
58 double Xbeta,
59 const Matrix<>& gamma){
60
61 Matrix<> cat_prob(1, ncat-1);
62 Matrix<> prob(1, ncat);
63 for (int j=0; j< ncat-1; ++j){
64 cat_prob(0, j) = pnorm(gamma[j+1] - Xbeta, 0.0, 1.0);
65 }
66 prob(0, ncat-1) = 1 - cat_prob(0, ncat-2);
67 prob(0, 0) = cat_prob(0, 0);
68 for (int j=1; j<(ncat-1); ++j){
69 prob(0, j) = cat_prob(0,j) - cat_prob(0, j-1);
70 }
71 double like = prob(0,Y-1);
72
73 return like;
74 }
75
oprobit_log_postLX(int j,const int ncat,const Matrix<> & gamma_p,const Matrix<> & gamma,const Matrix<> & Y,const Matrix<> & X,const Matrix<> & beta,const Matrix<> & tune,const int gammafixed)76 static double oprobit_log_postLX( int j,
77 const int ncat,
78 const Matrix<>& gamma_p,
79 const Matrix<>& gamma,
80 const Matrix<>& Y,
81 const Matrix<>& X,
82 const Matrix<>& beta,
83 const Matrix<>& tune,
84 const int gammafixed){
85 const int N = Y.rows();
86 double loglikerat = 0.0;
87 double loggendenrat = 0.0;
88 Matrix<> Xbeta = X*t(beta(j,_));
89
90 if (gammafixed==1){
91 for ( int i=0; i<N; ++i){
92 int yi = Y[i];
93 if (yi == ncat){
94 loglikerat = loglikerat
95 + log(1.0 - pnorm(gamma_p(yi-1) - Xbeta[i], 0.0, 1.0) )
96 - log(1.0 - pnorm(gamma(yi-1) - Xbeta[i], 0.0, 1.0) );
97 }
98 else if (yi == 1){
99 loglikerat = loglikerat + log(pnorm(gamma_p(yi) - Xbeta[i], 0.0, 1.0) )
100 - log(pnorm(gamma(yi) - Xbeta[i], 0.0, 1.0) );
101 }
102 else{
103 loglikerat = loglikerat
104 + log(pnorm(gamma_p(yi) - Xbeta[i], 0.0, 1.0) -
105 pnorm(gamma_p(yi-1) - Xbeta[i], 0.0, 1.0) )
106 - log(pnorm(gamma(yi) - Xbeta[i], 0.0, 1.0) -
107 pnorm(gamma(yi - 1) - Xbeta[i], 0.0, 1.0) );
108 }
109 }
110 for ( int k =2; k <ncat; ++k){
111 loggendenrat = loggendenrat +
112 log(pnorm(gamma(k+1), gamma(k), tune[j]) -
113 pnorm(gamma_p(k-1), gamma(k), tune[j])) -
114 log(pnorm(gamma_p(k+1), gamma_p(k), tune[j]) +
115 pnorm(gamma(k-1), gamma_p(k), tune[j]));
116 }
117 }
118 else{
119 for ( int i=0; i<N; ++i){
120 int yi = Y[i];
121 if (yi == ncat){
122 loglikerat = loglikerat
123 + log(1.0 - pnorm(gamma_p(j, yi-1) - Xbeta[i], 0.0, 1.0) )
124 - log(1.0 - pnorm(gamma(j, yi-1) - Xbeta[i], 0.0, 1.0) );
125 }
126 else if (yi == 1){
127 loglikerat = loglikerat + log(pnorm(gamma_p(j, yi) - Xbeta[i], 0.0, 1.0) )
128 - log(pnorm(gamma(j, yi) - Xbeta[i], 0.0, 1.0) );
129 }
130 else{
131 loglikerat = loglikerat
132 + log(pnorm(gamma_p(j, yi) - Xbeta[i], 0.0, 1.0) -
133 pnorm(gamma_p(j, yi-1) - Xbeta[i], 0.0, 1.0) )
134 - log(pnorm(gamma(j, yi) - Xbeta[i], 0.0, 1.0) -
135 pnorm(gamma(j, yi - 1) - Xbeta[i], 0.0, 1.0) );
136 }
137 }
138 for ( int k =2; k <ncat; ++k){
139 loggendenrat = loggendenrat +
140 log(pnorm(gamma(j ,k+1), gamma(j, k), tune[j]) -
141 pnorm(gamma_p(j, k-1), gamma(j, k), tune[j])) -
142 log(pnorm(gamma_p(j, k+1), gamma_p(j, k), tune[j]) +
143 pnorm(gamma(j, k-1), gamma_p(j, k), tune[j]));
144 }
145 }
146 return loglikerat + loggendenrat;
147 }
148
149 template <typename RNGTYPE>
gaussian_ordinal_state_sampler_fixedsigma(rng<RNGTYPE> & stream,const int m,const Matrix<> & Y,const Matrix<> & X,const Matrix<> & beta,const Matrix<> & Sigma,const Matrix<> & P)150 Matrix<> gaussian_ordinal_state_sampler_fixedsigma(rng<RNGTYPE>& stream,
151 const int m,
152 const Matrix<>& Y,
153 const Matrix<>& X,
154 const Matrix<>& beta,
155 const Matrix<>& Sigma,
156 const Matrix<>& P){
157
158 const int ns = m + 1;
159 const int n = Y.rows();
160 Matrix<> F(n, ns);
161 Matrix<> pr1(ns, 1);
162 pr1[0] = 1;
163 Matrix<> py(ns, 1);
164 Matrix<> pstyt1(ns, 1);
165
166 for (int t=0; t<n ; ++t){
167 Matrix<> mu = X(t,_)*::t(beta);
168 for (int j = 0; j< ns; ++j){
169 py[j] = dnorm(Y[t], mu[j], sqrt(Sigma[0]));
170 }
171 if (t==0) pstyt1 = pr1;
172 else {
173 pstyt1 = ::t(F(t-1,_)*P);
174 }
175 Matrix<> unnorm_pstyt = pstyt1%py;
176 const Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
177 for (int j=0; j<ns ; ++j) F(t,j) = pstyt(j);
178
179 }
180
181 Matrix<int> s(n, 1);
182 Matrix<> ps = Matrix<>(n, ns);
183 ps(n-1,_) = F(n-1,_);
184 s(n-1) = ns;
185
186 Matrix<> pstyn = Matrix<>(ns, 1);
187 double pone = 0.0;
188 int t = n-2;
189 while (t >= 0){
190 int st = s(t+1);
191 Matrix<> Pst_1 = ::t(P(_,st-1));
192 Matrix<> unnorm_pstyn = F(t,_)%Pst_1;
193 pstyn = unnorm_pstyn/sum(unnorm_pstyn);
194 if (st==1) s(t) = 1;
195 else{
196 pone = pstyn(st-2);
197 if(stream.runif() < pone) s(t) = st-1;
198 else s(t) = st;
199 }
200 ps(t,_) = pstyn;
201 --t;
202 }
203 Matrix<> Sout(n, ns+1);
204 Sout(_, 0) = s(_,0);
205 for (int j = 0; j<ns; ++j){
206 Sout(_,j+1) = ps(_, j);
207 }
208 return Sout;
209 }
210
211
212 ////////////////////////////////////////////
213 // cMCMCoprobitChangeLinearAppxpoint implementation.
214 ////////////////////////////////////////////
215 template <typename RNGTYPE>
MCMCoprobitChange_impl(rng<RNGTYPE> & stream,const int m,const int ncat,const Matrix<> & Y,const Matrix<> & X,Matrix<> & beta,Matrix<> & beta_linear,Matrix<> & gamma,Matrix<> & P,Matrix<> & Sigma,Matrix<> & b0,Matrix<> & B0,const Matrix<> & A0,int burnin,int mcmc,int thin,int verbose,const Matrix<> & tune,bool chib,bool gammafixed,Matrix<> & beta_store,Matrix<> & beta_linear_store,Matrix<> & gamma_store,Matrix<> & Z_store,Matrix<> & P_store,Matrix<> & ps_store,Matrix<int> & s_store,double & logmarglike,double & loglike)216 void MCMCoprobitChange_impl(rng<RNGTYPE>& stream,
217 const int m, const int ncat,
218 const Matrix<>& Y, const Matrix<>& X,
219 Matrix<>& beta, Matrix<>& beta_linear, Matrix<>& gamma, Matrix<>& P,
220 Matrix<>& Sigma,
221 Matrix<>& b0, Matrix<>& B0, const Matrix<>& A0,
222 int burnin, int mcmc, int thin,
223 int verbose, const Matrix<>& tune, // const Matrix<int>& tdf,
224 bool chib, bool gammafixed,
225 Matrix<>& beta_store, Matrix<>& beta_linear_store,
226 Matrix<>& gamma_store, Matrix<>& Z_store,
227 Matrix<>& P_store, Matrix<>& ps_store, Matrix<int>& s_store,
228 double& logmarglike, double& loglike)
229 {
230 const int tot_iter = burnin + mcmc;
231 const int nstore = mcmc / thin;
232 const int N = Y.rows();
233 const int ns = m + 1;
234 const int k = X.cols();
235 const int gk = ncat + 1;
236 const Matrix<> B0inv = invpd(B0);
237 Matrix<> Z(N, 1);
238 Matrix<> accepts(ns, 1);
239 Matrix<> gamma_p = gamma;
240
241 //MCMC loop
242 int count = 0;
243 for (int iter = 0; iter < tot_iter; ++iter){
244 // 1. Sample s
245 Matrix<> Sout = gaussian_ordinal_state_sampler_fixedsigma(stream, m, Y, X, beta_linear, Sigma, P);
246 Matrix<int> s = Sout(_, 0);
247 Matrix <double> ps(N, ns);
248 for (int j = 0; j<ns; ++j){
249 ps(_,j) = Sout(_,j+1);
250 }
251 // 2. Sample Z
252 for ( int i = 0; i<N; ++i) {
253 Matrix<> mu = X(i,_)*t(beta(s[i]-1,_));
254 int yi = Y[i];
255 Z[i] = stream.rtnorm_combo(mu[0], 1.0, gamma(s[i]-1, yi-1), gamma(s[i]-1, yi));
256 }
257
258 // 3. Sample beta
259 int beta_count = 0;
260 Matrix<int> beta_count_storage(ns, 1);
261 Matrix<int> nstate(ns, 1);
262 for (int j = 0; j<ns; ++j){
263 for (int i = 0; i<N; ++i){
264 if (s[i] == j + 1) {
265 nstate[j] = nstate[j] + 1;
266 }
267 }
268 beta_count = beta_count + nstate[j];
269 Matrix<> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
270 Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
271 Matrix<> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
272 Matrix<> XpX = t(Xj)*Xj;
273 Matrix<> XpZ = t(Xj)*Zj;
274 Matrix<> XpY = t(Xj)*yj;
275 Matrix<> Bn = invpd(B0 + XpX/Sigma[0]);
276 Matrix<> bn = Bn*(B0*b0 + XpY/Sigma[0]);
277 beta_linear(j,_) = stream.rmvnorm(bn, Bn);
278 Matrix<> Bn2 = invpd(B0 + XpX);
279 Matrix<> bn2 = Bn2*(B0*b0 + XpZ);
280 beta(j,_) = stream.rmvnorm(bn2, Bn2);
281 beta_count_storage[j] = beta_count;
282 }
283
284 // 4. Sample gamma
285 for (int j = 0; j < ns ; ++j){
286 for (int i = 2; i< ncat; ++i){
287 if (i==(ncat-1)){
288 gamma_p(j, i) = stream.rtbnorm_combo(gamma(j, i), ::pow(tune[j], 2.0),
289 gamma_p(j, i-1));
290 }
291 else {
292 gamma_p(j, i) = stream.rtnorm_combo(gamma(j, i), ::pow(tune[j], 2.0),
293 gamma_p(j, i-1), gamma(j, i+1));
294 }
295 }
296 Matrix<> Yj = Y((beta_count_storage[j] - nstate[j]), 0, (beta_count_storage[j] - 1), 0);
297 Matrix<> Xj = X((beta_count_storage[j] - nstate[j]), 0, (beta_count_storage[j] - 1), k-1);
298 double alpha = oprobit_log_postLX(j, ncat, gamma_p, gamma, Yj, Xj, beta, tune, gammafixed);
299
300 if (stream.runif() <= exp(alpha)){
301 gamma(j,_) = gamma_p(j,_);
302 accepts[j] = accepts[j] + 1;
303 }
304 }
305
306 // 5. Sample P
307 double shape1 = 0;
308 double shape2 = 0;
309 P(ns-1, ns-1) = 1;
310
311 for (int j =0; j< (ns-1); ++j){
312 shape1 = A0(j,j) + (double)nstate[j] - 1;
313 shape2 = A0(j,j+1) + 1; //
314 P(j,j) = stream.rbeta(shape1, shape2);
315 P(j,j+1) = 1 - P(j,j);
316 }
317
318 // load draws into sample array
319 if (iter >= burnin && ((iter % thin)==0)){
320 Matrix<> tbeta = ::t(beta);
321 Matrix<> tbetaLX = ::t(beta_linear);
322 for (int i=0; i<(ns*k); ++i){
323 beta_store(count,i) = tbeta[i];
324 beta_linear_store(count,i) = tbetaLX[i];
325 }
326 Matrix<> tgamma = ::t(gamma);
327 for (int i=0; i<(ns*gk); ++i)
328 gamma_store(count,i) = tgamma[i];
329 for (int j=0; j<ns*ns; ++j)
330 P_store(count,j)= P[j];
331 for (int l=0; l<N ; ++l)
332 ps_store(l,_) = ps_store(l,_) + ps(l,_);
333 s_store(count,_) = s;
334 Z_store(count,_) = Z;
335
336 ++count;
337
338 } // end of if(iter...)
339
340
341 // print output to stdout
342 if(iter > 1 && verbose > 0 && iter % verbose == 0){
343 Rprintf("\n\nMCMCoprobitChange iteration %i of %i \n", (iter+1), tot_iter);
344 for (int j = 0;j<ns; ++j){
345 double acceptancerate = accepts[j]/iter;
346 Rprintf("\n\n Acceptance rate for state %i is %10.5f \n", j+1, acceptancerate);
347 }
348 for (int j = 0;j<ns; ++j){
349 Rprintf("\n The number of observations in state %i is %10.5d", j+1, nstate[j]);
350 }
351 for (int j = 0; j<ns; ++j){
352 Rprintf("\n beta %i = ", j);
353 for (int i = 0; i<k; ++i){
354 Rprintf("%10.5f", beta(j, i));
355 }
356 }
357 // for (int j = 0; j<ns; ++j){
358 // Rprintf("\n beta_linear %i = ", j);
359 // for (int i = 0; i<k; ++i){
360 // Rprintf("%10.5f", beta_linear(j, i));
361 // }
362 // }
363 for (int j = 0; j<ns; ++j){
364 Rprintf("\n gamma %i = ", j);
365 for (int i = 0; i<(ncat-2); ++i){
366 Rprintf("%10.5f", gamma(j, i+2));
367 }
368 }
369 }
370
371 R_CheckUserInterrupt();
372
373 }// end MCMC loop
374
375 if(chib==1){
376 Matrix<> betast = meanc(beta_store);
377 Matrix<> betastLX = meanc(beta_linear_store);
378 Matrix<double, Row> beta_st(ns, k);
379 Matrix<double, Row> beta_linear_st(ns, k);
380 for (int j = 0; j<ns*k; ++j){
381 beta_st[j] = betast[j];
382 beta_linear_st[j] = betastLX[j];
383 }
384 Matrix<> gammast = meanc(gamma_store);
385 Matrix<double, Row> gamma_st(ns, gk);
386 for (int j = 0; j<ns*gk; ++j){
387 gamma_st[j] = gammast[j];
388 }
389
390 // for (int j = 0; j<ns; ++j){
391 // Rprintf("\n gamma_st %i = ", j);
392 // for (int i = 0; i<gk; ++i){
393 // Rprintf("%10.5f", gamma_st(j, i));
394 // }
395 // }
396
397 Matrix<> P_vec_st = meanc(P_store);
398 const Matrix<> P_st(ns, ns);
399 for (int j = 0; j< ns*ns; ++j){
400 P_st[j] = P_vec_st[j];
401 }
402 // storage
403 Matrix<> pdf_numer_store(nstore, 1);
404 Matrix<> pdf_alpha_store(nstore, 1);
405 Matrix<> pdf_P_store(nstore, ns);
406
407 // 1. gamma
408 Matrix<> densityq(nstore, ns);
409 Matrix<> alpha(nstore, ns);
410 for (int iter = 0; iter < nstore; ++iter){
411 int beta_count = 0;
412 Matrix<int> nstate(ns, 1);
413
414 Matrix<double, Row> gamma_g(ns, gk);
415 for (int h = 0; h<(ns*gk); ++h){
416 gamma_g[h] = gamma_store(iter, h);
417 }
418
419 Matrix<double, Row> beta_g(ns, k);
420 for (int h = 0; h<(ns*k); ++h){
421 beta_g[h] = beta_store(iter, h);
422 }
423 Matrix<> pdf_numer(ns, 1);
424 for (int j = 0; j <ns ; ++j){
425 for (int i = 0; i<N; ++i){
426 if (s_store(iter, i) == (j+1)) {
427 nstate[j] = nstate[j] + 1;
428 }
429 }
430 beta_count = beta_count + nstate[j];
431
432 Matrix<int> Yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
433 Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
434 pdf_numer(j) = oprobit_log_postLX(j, ncat, gamma_st, gamma_g, Yj, Xj, beta_g,
435 tune, gammafixed);
436 for (int h = 2; h<ncat; ++h){
437 if (h == (ncat-1)){
438 densityq(iter, j) = densityq(iter, j) + log(dtnormLX(gamma_st(j, h), gamma_g(j, h), tune[j],
439 gamma_st(j, h-1), 300));
440 }
441 else {
442 densityq(iter, j) = densityq(iter, j) + log(dtnormLX(gamma_st(j, h), gamma_g(j, h), tune[j],
443 gamma_st(j, h-1), gamma_g(j, h+1)));
444 }
445 }
446 }
447
448 if (sum(pdf_numer) > 0){
449 pdf_numer_store(iter) = 0;
450 }
451 else{
452 pdf_numer_store(iter) = sum(pdf_numer);
453 }
454
455 }
456 double numerator = sum(meanc(pdf_numer_store)) + sum(meanc(densityq));
457
458 for (int iter = 0; iter < nstore; ++iter){
459 Matrix<> Sout = gaussian_ordinal_state_sampler_fixedsigma(stream, m, Y, X, beta_linear, Sigma, P);
460 Matrix<int> s = Sout(_, 0);
461 for ( int i = 0; i<N; ++i) {
462 const Matrix<> mu = X(i,_)*t(beta(s[i]-1,_));
463 Z[i] = stream.rtnorm_combo(mu[0], 1.0, gamma_st(s[i]-1, Y[i]-1), gamma_st(s[i]-1, Y[i]));
464 }
465 int beta_count = 0;
466 Matrix<int> beta_count_storage(ns, 1);
467 Matrix<int> nstate(ns, 1);
468 for (int j = 0; j <ns ; ++j){
469 for (int i = 0; i<N; ++i){
470 if (s[i] == (j+1)) {
471 nstate[j] = nstate[j] + 1;
472 }
473 }
474 beta_count = beta_count + nstate[j];
475 Matrix<> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
476 Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
477 Matrix<> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
478 Matrix<> XpX = t(Xj)*Xj;
479 Matrix<> XpZ = t(Xj)*Zj;
480 Matrix<> XpY = t(Xj)*yj;
481 Matrix<> Bn = invpd(B0 + XpX/Sigma[0]);
482 Matrix<> bn = Bn*(B0*b0 + XpY/Sigma[0]);
483 beta_linear(j,_) = stream.rmvnorm(bn, Bn);
484 Matrix<> Bn2 = invpd(B0 + XpX);
485 Matrix<> bn2 = Bn2*(B0*b0 + XpZ);
486 beta(j,_) = stream.rmvnorm(bn2, Bn2);
487 beta_count_storage[j] = beta_count;
488 }
489 // Sample P
490 double shape1 = 0;
491 double shape2 = 0;
492 P(ns-1, ns-1) = 1;
493
494 for (int j =0; j< (ns-1); ++j){
495 shape1 = A0(j,j) + nstate[j] - 1;
496 shape2 = A0(j,j+1) + 1; //
497 P(j,j) = stream.rbeta(shape1, shape2);
498 P(j,j+1) = 1 - P(j,j);
499 }
500 Matrix<> alpha(ns, 1);
501 for (int j = 0; j < ns ; ++j){
502 for (int i = 2; i< ncat; ++i){
503 if (i==(ncat-1)){
504 gamma_p(j, i) = stream.rtbnorm_combo(gamma_st(j, i), ::pow(tune[j], 2.0),
505 gamma_p(j, i-1));
506 }
507 else {
508 gamma_p(j, i) = stream.rtnorm_combo(gamma_st(j, i), ::pow(tune[j], 2.0),
509 gamma_p(j, i-1), gamma_st(j, i+1));
510 }
511 }
512 Matrix<int> Yj = Y((beta_count_storage[j] - nstate[j]), 0, (beta_count_storage[j] - 1), 0);
513 Matrix<> Xj = X((beta_count_storage[j] - nstate[j]), 0, (beta_count_storage[j] - 1), k-1);
514 alpha[j] = oprobit_log_postLX(j, ncat, gamma_p, gamma_st, Yj, Xj, beta, tune, gammafixed);
515 }
516 if (sum(alpha) > 0){
517 pdf_alpha_store(iter) = 0;
518 }
519 else{
520 pdf_alpha_store(iter) = sum(alpha);
521 }
522
523 }
524 double denominator = mean(pdf_alpha_store);
525 double pdf_gamma = numerator - denominator;
526
527 // 2. beta
528 Matrix<> density_beta(nstore, ns);
529 for (int iter = 0; iter < nstore; ++iter){
530 Matrix<> Sout = gaussian_ordinal_state_sampler_fixedsigma(stream, m, Y, X, beta_linear, Sigma, P);
531 Matrix<int> s = Sout(_, 0);
532 for ( int i = 0; i<N; ++i) {
533 const Matrix<> mu = X(i,_)*t(beta(s[i]-1,_));
534 Z[i] = stream.rtnorm_combo(mu[0], 1.0, gamma_st(s[i]-1, Y[i]-1), gamma_st(s[i]-1, Y[i]));
535 }
536 int beta_count = 0;
537 Matrix<int> beta_count_storage(ns, 1);
538 Matrix<int> nstate(ns, 1);
539 for (int j = 0; j <ns ; ++j){
540 for (int i = 0; i<N; ++i){
541 if (s[i] == (j+1)) {
542 nstate[j] = nstate[j] + 1;
543 }
544 }
545 beta_count = beta_count + nstate[j];
546 Matrix<> yj = Y((beta_count - nstate[j]), 0, (beta_count - 1), 0);
547 Matrix<> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
548 Matrix<> Zj = Z((beta_count - nstate[j]), 0, (beta_count - 1), 0);
549 Matrix<> XpX = t(Xj)*Xj;
550 Matrix<> XpZ = t(Xj)*Zj;
551 Matrix<> XpY = t(Xj)*yj;
552 Matrix<> Bn = invpd(B0 + XpX/Sigma[0]);
553 Matrix<> bn = Bn*(B0*b0 + XpY/Sigma[0]);
554 beta_linear(j,_) = stream.rmvnorm(bn, Bn);
555 Matrix<> Bn2 = invpd(B0 + XpX);
556 Matrix<> bn2 = Bn2*(B0*b0 + XpZ);
557 beta(j,_) = stream.rmvnorm(bn2, Bn2);
558 beta_count_storage[j] = beta_count;
559 density_beta(iter, j) = exp(lndmvn(t(beta_st(j,_)), bn2, Bn2));
560 }
561
562 // Sample P
563 double shape1 = 0;
564 double shape2 = 0;
565 P(ns-1, ns-1) = 1;
566
567 for (int j =0; j< (ns-1); ++j){
568 shape1 = A0(j,j) + nstate[j] - 1;
569 shape2 = A0(j,j+1) + 1; //
570 P(j,j) = stream.rbeta(shape1, shape2);
571 P(j,j+1) = 1 - P(j,j);
572 }
573 }
574
575 double pdf_beta = log(prod(meanc(density_beta)));
576
577 // 3. P
578 Matrix<> density_P(nstore, ns);
579 for (int iter = 0; iter < nstore; ++iter){
580 Matrix<> Sout = gaussian_ordinal_state_sampler_fixedsigma(stream, m, Y, X, beta_linear_st, Sigma, P);
581 Matrix <double> s = Sout(_, 0);
582 double shape1 = 0;
583 double shape2 = 0;
584 P(ns-1, ns-1) = 1;
585 Matrix <double> P_addN(ns, 1);
586 for (int j = 0; j<ns; ++j){
587 for (int i = 0; i<N; ++i){
588 if (s[i] == (j+1)) {
589 P_addN[j] = P_addN[j] + 1;
590 }
591 }
592 }
593
594 for (int j =0; j< (ns-1); ++j){
595 shape1 = A0(j,j) + P_addN[j] - 1;
596 shape2 = A0(j,j+1) + 1;
597 P(j,j) = stream.rbeta(shape1, shape2);
598 P(j,j+1) = 1 - P(j,j);
599 density_P(iter, j) = dbeta(P_st(j,j), shape1, shape2);
600 }
601 density_P(iter, ns-1) = 1;
602
603 }
604 double pdf_P = log(prod(meanc(density_P)));
605
606 // likelihood
607 Matrix<> F = Matrix<>(N, ns);
608 Matrix<> like(N, 1);
609 Matrix<> pr1 = Matrix<>(ns, 1);
610 pr1[0] = 1;
611 Matrix<> py(ns, 1);
612 Matrix<> pstyt1(ns, 1);
613
614 for (int t=0; t< N ; ++t){
615 Matrix<> mu = X(t,_)*::t(beta_st);
616 for (int j=0; j<ns ; ++j){
617 int yt = Y[t];
618 py[j] = oprobit_pdfLX(ncat, yt, mu[j], gamma_st(j,_));
619 }
620 if (t==0) pstyt1 = pr1;
621 else {
622 pstyt1 = ::t(F(t-1,_)*P_st);
623 }
624 Matrix<> unnorm_pstyt = pstyt1%py;
625 Matrix<> pstyt = unnorm_pstyt/sum(unnorm_pstyt);
626 for (int j=0; j<ns ; ++j){
627 F(t,j) = pstyt(j);
628 }
629 like[t] = sum(unnorm_pstyt);
630 }
631
632 loglike = sum(log(like));
633
634 // log prior ordinate
635 Matrix<> density_beta_prior(ns, 1);
636 Matrix<> density_P_prior(ns, 1);
637 density_P[ns-1] = 1; //
638
639 for (int j=0; j<ns ; ++j){
640 density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
641 }
642
643 for (int j =0; j< (ns-1); ++j){
644 density_P_prior[j] = log(dbeta(P_st(j,j), A0(j,j), A0(j,j+1)));
645 }
646 double density_gamma_prior = ns*(ncat-2)*log(dunif(1, 0, 10));
647
648 double logprior = sum(density_beta_prior) + sum(density_P_prior) + density_gamma_prior;
649 logmarglike = (loglike + logprior) - (pdf_beta + pdf_P + pdf_gamma);
650
651 if (verbose >0 ){
652 Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
653 Rprintf("loglike = %10.5f\n", loglike);
654 Rprintf("log_prior = %10.5f\n", logprior);
655 Rprintf("log_beta = %10.5f\n", pdf_beta);
656 Rprintf("log_P = %10.5f\n", pdf_P);
657 Rprintf("log_gamma = %10.5f\n", pdf_gamma);
658 }
659 } // end of marginal likelihood
660
661 }//end
662
663 extern "C"{
cMCMCoprobitChange(double * betaout,double * betalinearout,double * gammaout,double * Pout,double * psout,double * sout,const double * Ydata,const double * Xdata,const int * Xrow,const int * Xcol,const int * m,const int * ncat,const int * burnin,const int * mcmc,const int * thin,const int * verbose,const double * tunedata,const int * uselecuyer,const int * seedarray,const int * lecuyerstream,const double * betastart,const double * betalinearstart,const double * gammastart,const double * Pstart,const double * sigmastart,const double * a,const double * b,const double * b0data,const double * B0data,const double * A0data,double * logmarglikeholder,double * loglikeholder,const int * chib,const int * gammafixed)664 void cMCMCoprobitChange(double *betaout, double *betalinearout,
665 double *gammaout, double *Pout, double *psout, double *sout,
666 const double *Ydata,
667 const double *Xdata, const int *Xrow, const int *Xcol,
668 const int *m, const int *ncat,
669 const int *burnin, const int *mcmc, const int *thin, const int *verbose,
670 const double *tunedata, // const int *tdfdata,
671 const int *uselecuyer, const int *seedarray, const int *lecuyerstream,
672 const double *betastart, const double *betalinearstart,
673 const double *gammastart, const double *Pstart,
674 const double *sigmastart, const double *a, const double *b,
675 const double *b0data, const double *B0data,
676 const double *A0data,
677 double *logmarglikeholder, double *loglikeholder,
678 const int *chib, const int *gammafixed)
679 {
680
681 // pull together Matrix objects
682 const Matrix<> Y(*Xrow, 1, Ydata);
683 const Matrix<> X(*Xrow, *Xcol, Xdata);
684 const int nstore = *mcmc / *thin;
685 const int N = *Xrow;
686 const int k = *Xcol;
687 const int gk = *ncat + 1;
688 const int ns = *m + 1;
689
690 // generate starting values
691 Matrix<> beta(ns, k, betastart);
692 Matrix<> beta_linear(ns, k, betalinearstart);
693 Matrix<> Sigma(1, 1, sigmastart);
694 Matrix<> P(ns, ns, Pstart);
695 Matrix<> b0(k, 1, b0data);
696 Matrix<> B0(k, k, B0data);
697 Matrix<> tune(ns, 1, tunedata);
698 Matrix<> A0(ns, ns, A0data);
699 double logmarglike;
700 double loglike;
701
702 // storage matrices
703 Matrix<> beta_store(nstore, ns*k);
704 Matrix<> beta_linear_store(nstore, ns*k);
705 Matrix<> Z_store(nstore, N);
706 Matrix<> P_store(nstore, ns*ns);
707 Matrix<> ps_store(N, ns);
708 Matrix<int> s_store(nstore, N);
709
710 Matrix<> gamma(ns, gk, gammastart);
711 Matrix<> gamma_store(nstore, ns*gk);
712 MCMCPACK_PASSRNG2MODEL(MCMCoprobitChange_impl, *m, *ncat,
713 Y, X, beta, beta_linear, gamma, P, Sigma,
714 b0, B0, A0,
715 *burnin, *mcmc, *thin, *verbose, tune,
716 *chib, *gammafixed,
717 beta_store, beta_linear_store, gamma_store, Z_store,
718 P_store, ps_store, s_store,
719 logmarglike, loglike);
720
721 logmarglikeholder[0] = logmarglike;
722 loglikeholder[0] = loglike;
723
724 // return output
725 for (int i = 0; i<(nstore*ns*k); ++i){
726 betaout[i] = beta_store[i];
727 betalinearout[i] = beta_linear_store[i];
728 }
729 for (int i = 0; i<(nstore*ns*gk); ++i){
730 gammaout[i] = gamma_store[i];
731 }
732 for (int i = 0; i<(nstore*ns*ns); ++i){
733 Pout[i] = P_store[i];
734 }
735 for (int i = 0; i<(N*ns); ++i){
736 psout[i] = ps_store[i];
737 }
738 for (int i = 0; i<(nstore*N); ++i){
739 sout[i] = s_store[i];
740 }
741 }
742 }
743 #endif
744
745
746