1 //////////////////////////////////////////////////////////////////////////
2 // MCMCfcds.h is the header file for MCMCfcds.cc. It contains declarations
3 // for a number of functions that produce draws from full conditional
4 // distributions.
5 //
6 // Andrew D. Martin
7 // Dept. of Political Science
8 // Washington University in St. Louis
9 // admartin@wustl.edu
10 //
11 // Kevin M. Quinn
12 // Dept. of Government
13 // Harvard University
14 // kevin_quinn@harvard.edu
15 //
16 // This software is distributed under the terms of the GNU GENERAL
17 // PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
18 // file for more information.
19 //
20 // KQ 6/10/2004
21 // DBP 7/01/2007 [ported to scythe 1.0.x (partial)]
22 //
23 // Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
24 // Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
25 // and Jong Hee Park
26 //////////////////////////////////////////////////////////////////////////
27
28
29
30 #ifndef MCMCFCDS_H
31 #define MCMCFCDS_H
32
33 #include "matrix.h"
34 #include "rng.h"
35 #include "stat.h"
36 #include "smath.h"
37 #include "ide.h"
38 #include "la.h"
39 #include "distributions.h"
40
41 #include <iostream>
42
43 using namespace std;
44 using namespace scythe;
45
46 // linear regression with Gaussian errors beta draw
47 // (multivariate Normal prior)
48 // regression model is y = X * beta + epsilon, epsilon ~ N(0,sigma2)
49 // XpX is X'X
50 // XpY is X'y
51 // b0 is the prior mean of beta
52 // B0 is the prior precision (the inverse variance) of beta
53 template <typename RNGTYPE>
54 Matrix<double>
NormNormregress_beta_draw(const Matrix<> & XpX,const Matrix<> & XpY,const Matrix<> & b0,const Matrix<> & B0,double sigma2,rng<RNGTYPE> & stream)55 NormNormregress_beta_draw (const Matrix<>& XpX, const Matrix<>& XpY,
56 const Matrix<>& b0, const Matrix<>& B0, double sigma2,
57 rng<RNGTYPE>& stream)
58 {
59 // this function gets the cross-product matrix X'X and the matrix X'Y
60 // to minimize the amount of computation within the function
61 const unsigned int k = XpX.cols ();
62 const double sig2_inv = 1.0 / sigma2;
63 const Matrix<> sig_beta = invpd (B0 + XpX * sig2_inv);
64 const Matrix<> C = cholesky (sig_beta);
65 const Matrix<> betahat = sig_beta * gaxpy(B0, b0, XpY*sig2_inv);
66
67 return( gaxpy(C, stream.rnorm(k,1, 0, 1), betahat) );
68 }
69
70 // linear regression with Gaussian errors sigma2 draw
71 // (inverse-Gamma prior)
72 // regression model is y = X * beta + epsilon, epsilon ~ N(0,sigma2)
73 // c0/2 is the prior shape parameter for sigma2
74 // d0/2 is the prior scale parameter for sigma2
75 template <typename RNGTYPE>
76 double
NormIGregress_sigma2_draw(const Matrix<> & X,const Matrix<> & Y,const Matrix<> & beta,double c0,double d0,rng<RNGTYPE> & stream)77 NormIGregress_sigma2_draw (const Matrix <> &X, const Matrix <> &Y,
78 const Matrix <> &beta, double c0, double d0,
79 rng<RNGTYPE>& stream)
80 {
81
82 const Matrix <> e = gaxpy(X, (-1*beta), Y);
83 const Matrix <> SSE = crossprod (e);
84 const double c_post = (c0 + X.rows ()) * 0.5;
85 const double d_post = (d0 + SSE[0]) * 0.5;
86
87 return stream.rigamma (c_post, d_post);
88 }
89
90 // Bayesian quantile (including median) regression beta draw
91 // (multivariate Normal prior)
92 //
93 // b0 is the prior mean of beta(tau)
94 // B0 is the prior precision (the inverse variance) of beta(tau)
95 template <typename RNGTYPE>
96 Matrix<double>
ALaplaceNormregress_beta_draw(double tau,const Matrix<> & X,const Matrix<> & Y,const Matrix<> & weights,const Matrix<> & b0,const Matrix<> & B0,rng<RNGTYPE> & stream)97 ALaplaceNormregress_beta_draw (double tau, const Matrix<>& X, const Matrix<>& Y, const Matrix<>& weights,
98 const Matrix<>& b0, const Matrix<>& B0,
99 rng<RNGTYPE>& stream)
100 {
101
102 const unsigned int k = X.cols();
103 const unsigned int n_obs = X.rows();
104 Matrix<> U(Y);
105 if (tau!=0.5){
106 U -= (1.0-2.0*tau)*weights;
107 }
108 Matrix<> XtwX(k,k,false);
109 Matrix<> XtwU(k,1,false);
110 double temp_x = 0.0;
111 double temp_u = 0.0;
112
113 //Calculate XtwU where w denotes a diagonal matrix with the augmented data (weights) on the diagonal
114 for (unsigned int i=0; i<k; ++i){
115 for (unsigned int m=0; m<n_obs; ++m){
116 temp_u += X(m,i)*U(m)/weights(m);
117 }
118 XtwU(i) = temp_u;
119 temp_u = 0.0;
120 }
121
122 //Calculate XtwX
123 for (unsigned int i=0; i<k; ++i){
124 for (unsigned int j=0; j<(i+1); ++j){
125 for (unsigned int m=0; m<n_obs; ++m){
126 temp_x += X(m,i)*X(m,j)/weights(m);
127 }
128 XtwX(i,j) = temp_x;
129 XtwX(j,i) = temp_x;
130 temp_x = 0.0;
131 }
132 }
133
134 const Matrix<> var_matrix_beta = invpd(B0+0.5*XtwX);
135 const Matrix<> C = cholesky(var_matrix_beta);
136 const Matrix<> betahat = var_matrix_beta*gaxpy(B0,b0,0.5*XtwU);
137
138 return( gaxpy(C, stream.rnorm(k,1, 0, 1), betahat) );
139 }
140
141 // This function draws from the full conditional distribution of the latent random variables (weights) under quantile regression (including median regression) and returns a column vector of those weights.
142
143 template <typename RNGTYPE>
144 Matrix<double>
ALaplaceIGaussregress_weights_draw(const Matrix<> & abse,rng<RNGTYPE> & stream)145 ALaplaceIGaussregress_weights_draw (const Matrix <> &abse,
146 rng<RNGTYPE>& stream)
147
148 {
149 const Matrix<double> nu_params = pow(abse,-1.0);
150 Matrix<> w(abse);
151 const unsigned int n_obs = abse.rows();
152
153 // The inverse Gaussian distribution
154
155 for (unsigned int i=0; i<n_obs; ++i){
156 double chisq = stream.rchisq(1);
157 double nu = nu_params(i);
158 double smallroot = nu*(nu*chisq+1.0-std::sqrt(nu*nu*chisq*chisq+2.0*nu*chisq));
159 unsigned int q = stream.rbern(nu/(nu+smallroot));
160 if (q == 1){
161 w(i) = smallroot;
162 }
163 else{
164 w(i) = nu*nu/smallroot;
165 }
166 }
167 return(pow(w,-1.0));
168 }
169
170 /////////////////////////////////////////////////////
171
172 // Functions for the quantile regression stochastic search variable selection
173
174 struct COV_TRIAL{
175 Matrix<> Cnew;
176 bool newtrial;
177 double logdetminhalf;
178 };
179
180 // updating the indicator variable corresponding to whether a
181 // covariate is included in the model, given that it was
182 // previously absent
183 template <typename RNGTYPE>
184 COV_TRIAL
QR_SSVS_covariate_trials_draw_absent(const Matrix<> & C,const Matrix<> & X_gamma,const Matrix<> & U,const Matrix<> & newXcol,unsigned int row_index,const Matrix<> & weights,double pi0,double newlambda,double logolddetminhalf,rng<RNGTYPE> & stream)185 QR_SSVS_covariate_trials_draw_absent(const Matrix<>& C, const Matrix<>& X_gamma, const Matrix<>& U,
186 const Matrix<>& newXcol, unsigned int row_index, const Matrix<>& weights, double pi0, double newlambda, double logolddetminhalf, rng<RNGTYPE>& stream){
187 const unsigned int n_obs = U.rows();
188 const unsigned int k = C.rows();
189
190 //Calculate new row required to update the Cholesky decomposition
191 Matrix<> XUXnewtXnew(k+1,1,false);
192 double temp_xux1 = 0.0;
193 double temp_xux2 = 0.0;
194 double temp_xux3 = 0.0;
195
196 //Calculate XUXnewtXnew
197 for (unsigned int i=0; i<k-1; ++i){
198 for (unsigned int m=0; m<n_obs; ++m){
199 temp_xux1 += X_gamma(m,i)*newXcol(m)/weights(m);
200 }
201 XUXnewtXnew(i) = 0.5*temp_xux1;
202 temp_xux1 = 0.0;
203 }
204 for (unsigned int m=0; m<n_obs; ++m){
205 temp_xux2 += U(m)*newXcol(m)/weights(m);
206 temp_xux3 += newXcol(m)*newXcol(m)/weights(m);
207 }
208 XUXnewtXnew(k-1) = 0.5*temp_xux2;
209 XUXnewtXnew(k) = 0.5*temp_xux3+newlambda;
210
211 //Obtain the Cholesky Decomposition of the new matrix formed by adding newXcol onto the right hand side
212
213 Matrix<> z(k,1,false);
214 for (unsigned int i = 0; i < k; ++i) {
215 double sum = 0;
216 for (unsigned int j = 0; j < i; ++j) {
217 sum += C(i,j) * z(j);
218 }
219 z(i) = (XUXnewtXnew(i) - sum) / C(i, i);
220 }
221 double rho = std::sqrt(XUXnewtXnew(k)-crossprod(z)(0));
222
223 Matrix<> Cnew(k+1, k+1, true, 0.0);
224 Cnew(0,0,k-1,k-1) = C;
225 Cnew(k,0,k,k-1) = z;
226 Cnew(k,k) = rho;
227
228 // Permuting the Cholesky decomposition so that it corresponds to what would be obtained if the X matrix included the covariate
229
230 Matrix<> temp(Cnew);
231 if (row_index != 0){
232 temp(0,0,row_index-1,k) = Cnew(0,0,row_index-1,k);
233 }
234 temp(row_index,_) = Cnew(k,_);
235 temp(row_index+1,0,k,k) = Cnew(row_index,0,k-1,k);
236
237 // Givens rotations
238
239 Matrix<> Q(2,2,false);
240 for (unsigned int i=k; i>row_index; --i)
241 {
242 double two_norm = std::sqrt(temp(row_index,i)*temp(row_index,i)
243 +temp(row_index,i-1)*temp(row_index,i-1));
244 Q(0,0) = temp(row_index,i-1)/two_norm;
245 Q(1,0) = temp(row_index,i)/two_norm;
246 Q(1,1) = Q(0,0);
247 Q(0,1) = -1.0*Q(1,0);
248 if (i!=k){
249 temp(i+1,i-1,k,i) = temp(i+1,i-1,k,i) * Q;
250 }
251 double temp2 = temp(i,i-1);
252 temp(i,i-1) = Q(0,0)*temp2;
253 temp(i,i) = Q(0,1)*temp2;
254 if (temp(i,i) < 0){
255 temp(i,i,k,i) = -1.0*temp(i,i,k,i);
256 }
257 temp(row_index,i-1) = two_norm;
258 temp(row_index,i) = 0.0;
259 }
260 Cnew=temp;
261
262 //Work out -0.5*log(det(Cnew'Cnew))
263 double lognewdetminhalf = 0.0;
264 for (unsigned int i=0; i<k; ++i){
265 lognewdetminhalf -= std::log(Cnew(i,i));
266 }
267
268
269
270 double log_g0 = logolddetminhalf-0.5*C(k-1,k-1)*C(k-1,k-1);
271 double log_g1 = 0.5*std::log(newlambda)+lognewdetminhalf-0.5*Cnew(k,k)*Cnew(k,k);
272
273 double log_ratio = log_g0+std::log(1.0-pi0)-log_g1-std::log(pi0);
274 double success_prob = 1.0/(1.0+std::exp(log_ratio));
275 bool new_covariate_trial = stream.rbern(success_prob);
276 COV_TRIAL result;
277 result.newtrial = new_covariate_trial;
278 if (new_covariate_trial == false){
279 result.Cnew = C;
280 result.logdetminhalf = logolddetminhalf;
281 }
282 else {
283 result.Cnew = Cnew;
284 result.logdetminhalf = lognewdetminhalf;
285 }
286 return result;
287 }
288
289 // updating the indicator variable corresponding to whether a
290 // covariate is included in the model, given that it was
291 // previously present in the model
292
293 template <typename RNGTYPE>
294 COV_TRIAL
QR_SSVS_covariate_trials_draw_present(const Matrix<> & C,unsigned int row_index,unsigned int n_obs,double pi0,double oldlambda,double logolddetminhalf,rng<RNGTYPE> & stream)295 QR_SSVS_covariate_trials_draw_present(const Matrix<>& C, unsigned int row_index, unsigned int n_obs, double pi0, double oldlambda, double logolddetminhalf, rng<RNGTYPE>& stream){
296 unsigned int k = C.rows();
297
298 // Permuting the Cholesky decomposition so that it corresponds to what would be obtained if the X matrix had the covariate in the final column
299
300 Matrix<> temp(C);
301 if (row_index != 0){
302 temp(0,0,row_index-1,k-1) = C(0,0,row_index-1,k-1);
303 }
304 temp(k-1,_) = C(row_index,_);
305 temp(row_index,0,k-2,k-1) = C(row_index+1,0,k-1,k-1);
306
307 // Givens rotations
308
309 Matrix<> Q(2,2,false);
310 for (unsigned int i=row_index; i<k-1; ++i)
311 {
312 double two_norm = std::sqrt(temp(i,i)*temp(i,i)
313 +temp(i,i+1)*temp(i,i+1));
314 Q(0,0) = temp(i,i)/two_norm;
315 Q(1,0) = temp(i,i+1)/two_norm;
316 Q(1,1) = Q(0,0);
317 Q(0,1) = -1.0*Q(1,0);
318 if (i!=k-2){
319 temp(i+1,i,k-2,i+1) = temp(i+1,i,k-2,i+1) * Q;
320 }
321 double temp2 = temp(k-1,i);
322 temp(k-1,i) = Q(0,0)*temp2;
323 temp(k-1,i+1) = Q(0,1)*temp2;
324 temp(i,i) = two_norm;
325 temp(i,i+1) = 0.0;
326 }
327 if (temp(k-1,k-1) < 0){
328 temp(k-1,k-1) = -1.0*temp(k-1,k-1);
329 }
330
331 Matrix<> Cnew = temp(0,0,k-2,k-2);
332
333 // Work out -1/2*log(det(Cnew'Cnew))
334 double lognewdetminhalf = 0.0;
335 for (unsigned int i=0; i<k-2; ++i){
336 lognewdetminhalf -= std::log(Cnew(i,i));
337 }
338
339
340 double log_g0 = lognewdetminhalf-0.5*Cnew(k-2,k-2)*Cnew(k-2,k-2);
341 double log_g1 = 0.5*std::log(oldlambda)+logolddetminhalf-0.5*C(k-1,k-1)*C(k-1,k-1);
342
343
344 double log_ratio = log_g0+std::log(1.0-pi0)-log_g1-std::log(pi0);
345 double success_prob = 1.0/(1.0+std::exp(log_ratio));
346 bool new_covariate_trial = stream.rbern(success_prob);
347 COV_TRIAL result;
348 result.newtrial = new_covariate_trial;
349 if (new_covariate_trial == false){
350 result.Cnew = Cnew;
351 result.logdetminhalf = lognewdetminhalf;
352 }
353 else {
354 result.Cnew = C;
355 result.logdetminhalf = logolddetminhalf;
356 }
357 return result;
358 }
359
360 // update betas using Cholesky decomposition
361 template <typename RNGTYPE>
362 Matrix<>
QR_SSVS_beta_draw(const Matrix<> & C,rng<RNGTYPE> & stream)363 QR_SSVS_beta_draw(const Matrix<>& C, rng<RNGTYPE>& stream){
364 unsigned int k = C.rows();
365 Matrix<> standnorm = stream.rnorm(k-1,1,0,1);
366 Matrix<> z(k-1,1,false);
367 z = t(C(k-1,0,k-1,k-2));
368 Matrix<> Q = z+standnorm*std::sqrt(2.0);
369 Matrix<> result(k-1,1,false);
370 for (int i = k-2; i >= 0; --i) {
371 double sum = 0;
372 for (unsigned int j = i+1; j < k-1; ++j) {
373 sum += C(j,i) * result(j);
374 }
375 result(i) = (Q(i) - sum) / C(i, i);
376 }
377 return result;
378 }
379
380 //hyperparameter pi0 updating
381
382 template <typename RNGTYPE>
383 double
QR_SSVS_pi0_draw(unsigned int n_uncert_cov,unsigned int tot_n_uncert_cov,double pi0a0,double pi0b0,rng<RNGTYPE> & stream)384 QR_SSVS_pi0_draw(unsigned int n_uncert_cov, unsigned int tot_n_uncert_cov,
385 double pi0a0, double pi0b0, rng<RNGTYPE>& stream){
386 double pi0a1 = pi0a0 + n_uncert_cov;
387 double pi0b1 = pi0b0 + tot_n_uncert_cov - n_uncert_cov;
388 return(stream.rbeta(pi0a1,pi0b1));
389 }
390
391 //update latent lambdas
392
393 template<typename RNGTYPE>
394 Matrix<double>
QR_SSVS_lambda_draw(const Matrix<> & beta_red,const Matrix<> & gamma,unsigned int tot_n_cov,unsigned int n_cert_cov,rng<RNGTYPE> & stream)395 QR_SSVS_lambda_draw(const Matrix<>& beta_red, const Matrix<>& gamma, unsigned int tot_n_cov, unsigned int n_cert_cov, rng<RNGTYPE>& stream)
396 {
397 unsigned int n_uncert_cov = tot_n_cov - n_cert_cov;
398 Matrix<> newlambda(n_uncert_cov,1,false);
399
400 for (unsigned int i=n_cert_cov; i<tot_n_cov; ++i){
401
402 unsigned int j = i-n_cert_cov;
403
404 if (gamma(i) == true){
405
406 unsigned int col_index = n_cert_cov;
407 //obtain column index of betas
408 for (unsigned int m=n_cert_cov; m<i; ++m){
409 if (gamma(m) == true){
410 ++col_index;
411 }
412 }
413
414 newlambda(j) = stream.rexp(0.5*(1.0+beta_red(col_index)*beta_red(col_index)));
415 }
416
417 else {
418 newlambda(j) = stream.rexp(0.5);
419 }
420
421 }
422 return newlambda;
423 }
424
425
426 /////////////////////////////////////////////////////
427
428 // update latent data for standard item response models
429 // only works for 1 dimensional case
430 template <typename RNGTYPE>
irt_Z_update1(Matrix<> & Z,const Matrix<int> & X,const Matrix<> & theta,const Matrix<> & eta,rng<RNGTYPE> & stream)431 void irt_Z_update1 (Matrix<>& Z, const Matrix<int>& X,
432 const Matrix<>& theta, const Matrix<>& eta, rng<RNGTYPE>& stream)
433 {
434 // define constants
435 const unsigned int J = theta.rows();
436 const unsigned int K = eta.rows();
437
438 // perform update from truncated Normal / standard Normals
439 for (unsigned int i = 0; i < J; ++i) {
440 for (unsigned int j = 0; j < K; ++j){
441 const double Z_mean = -eta(j,0) + theta(i) * eta(j,1);
442 if (X(i,j) == 1) {
443 Z(i,j) = stream.rtbnorm_combo(Z_mean, 1.0, 0);
444 } else if (X(i,j) == 0) {
445 Z(i,j) = stream.rtanorm_combo(Z_mean, 1.0, 0);
446 } else {
447 Z(i,j) = stream.rnorm(Z_mean, 1.0);
448 }
449 }
450 }
451 }
452
453 // update item (case, roll call) parameters for item response model
454 // note: works only for one-dimensional case
455 template <typename RNGTYPE>
irt_eta_update1(Matrix<> & eta,const Matrix<> & Z,const Matrix<> & theta,const Matrix<> & AB0,const Matrix<> & AB0ab0,rng<RNGTYPE> & stream)456 void irt_eta_update1 (Matrix<>& eta, const Matrix<>& Z,
457 const Matrix<>& theta, const Matrix<>& AB0, const Matrix<>& AB0ab0,
458 rng<RNGTYPE>& stream)
459 {
460
461 // define constants
462 const unsigned int J = theta.rows();
463 const unsigned int K = Z.cols();
464
465 // perform update
466 //const Matrix<double> Ttheta_star = t(cbind(-1.0*ones<double>(J,1),theta)); // only needed for option 2
467 const Matrix<> tpt(2,2);
468 for (unsigned int i = 0; i < J; ++i) {
469 const double theta_i = theta(i);
470 tpt(0,1) -= theta_i;
471 tpt(1,1) += std::pow(theta_i, 2.0);
472 }
473 tpt(1,0) = tpt(0,1);
474 tpt(0,0) = J;
475 const Matrix<> eta_post_var = invpd(tpt + AB0);
476 const Matrix<> eta_post_C = cholesky(eta_post_var);
477
478 for (unsigned int k = 0; k < K; ++k) {
479 const Matrix<> TZ(2, 1);
480 for (unsigned int j = 0; j < J; ++j) {
481 TZ[0] -= Z(j,k);
482 TZ[1] += Z(j,k) * theta[j];
483 }
484 const Matrix<> eta_post_mean = eta_post_var * (TZ + AB0ab0);
485 const Matrix<> new_eta = gaxpy(eta_post_C, stream.rnorm(2, 1, 0, 1),
486 eta_post_mean);
487 eta(k,0) = new_eta(0);
488 eta(k,1) = new_eta(1);
489 }
490 }
491
492 // update ability parameters (ideal points) for one dimensional
493 // item response model
494 // note: works only for one-dimensional case
495 template <typename RNGTYPE>
irt_theta_update1(Matrix<> & theta,const Matrix<> & Z,const Matrix<> & eta,double t0,double T0,const Matrix<> & theta_eq,const Matrix<> & theta_ineq,rng<RNGTYPE> & stream)496 void irt_theta_update1 (Matrix<>& theta, const Matrix<>& Z,
497 const Matrix<>& eta, double t0, double T0, const Matrix<>& theta_eq,
498 const Matrix<>& theta_ineq, rng<RNGTYPE>& stream)
499 {
500
501 const unsigned int J = Z.rows();
502 const unsigned int K = Z.cols();
503
504 // perform update from multivariate Normal
505 const double T0t0 = T0*t0;
506 const Matrix<> alpha = eta(_, 0);
507 const Matrix<> beta = eta(_, 1);
508 //const Matrix<double> tbeta = t(beta); // only neede for option 2
509 //const Matrix<double> talpha = t(alpha); // only needed for option 2
510
511 // calculate the posterior variance outside the justice specific loop
512 double theta_post_var = T0;
513 for (unsigned int i = 0; i < K; ++i){
514 theta_post_var += std::pow(beta(i), 2.0);
515 }
516 theta_post_var = 1.0 / theta_post_var;
517 const double theta_post_sd = std::sqrt(theta_post_var);
518
519 // sample for each justice
520 for (unsigned int j = 0; j < J; ++j) {
521 // no equality constraints
522 if (theta_eq(j) == -999) {
523 double betaTZjalpha = 0;
524 for (unsigned int k = 0; k < K; ++k){
525 betaTZjalpha += beta(k) * (Z(j,k) + alpha(k));
526 }
527 const double theta_post_mean = theta_post_var*(T0t0 + betaTZjalpha);
528 if (theta_ineq(j) == 0) { // no inequality constraint
529 theta(j) = theta_post_mean + stream.rnorm(0.0, theta_post_sd);
530 } else if (theta_ineq(j) > 0) { // theta[j] > 0
531 theta(j) = stream.rtbnorm_combo(theta_post_mean, theta_post_var, 0);
532 } else { // theta[j] < 0
533 theta(j) = stream.rtanorm_combo(theta_post_mean, theta_post_var, 0);
534 }
535 } else { // equality constraints
536 theta(j) = theta_eq(j);
537 }
538 }
539 }
540
541
542
543
544
545
546
547 // update latent Ystar values in paired comparisons model
548 template <typename RNGTYPE>
paircompare_Ystar_update(Matrix<> & Ystar,const Matrix<unsigned int> & MD,const Matrix<> & theta,const Matrix<> & alpha,rng<RNGTYPE> & stream)549 void paircompare_Ystar_update (Matrix<>& Ystar,
550 const Matrix<unsigned int>& MD,
551 const Matrix<>& theta, const Matrix<>& alpha,
552 rng<RNGTYPE>& stream){
553
554 const unsigned int N = MD.rows();
555 for (unsigned int i = 0; i < N; ++i){
556 const double alpha_i = alpha(MD(i,0));
557 const double theta_1 = theta(MD(i,1));
558 const double theta_2 = theta(MD(i,2));
559 const double mean = alpha_i * (theta_1 - theta_2);
560 const bool above = MD(i,1) == MD(i,3); // cand 1 chosen
561 const bool below = MD(i,2) == MD(i,3); // cand 2 chosen
562 if (above){
563 Ystar(i) = stream.rtbnorm_combo(mean, 1.0, 0.0);
564 }
565 else if (below){
566 Ystar(i) = stream.rtanorm_combo(mean, 1.0, 0.0);
567 }
568 else{
569 Ystar(i) = stream.rnorm(mean, 1.0);
570 }
571 }
572 }
573
574
575
576
577
578 // update theta values in paired comparisons model
579 template <typename RNGTYPE>
paircompare_theta_update(Matrix<> & theta,const Matrix<> & Ystar,const Matrix<unsigned int> & MD,const Matrix<> & alpha,const Matrix<unsigned int> & theta_n,const Matrix<> & theta_eq,const Matrix<> & theta_ineq,const vector<vector<double * >> & theta_Ystar_ptr,const vector<vector<double * >> & theta_alpha_ptr,const vector<vector<double * >> & theta_theta_ptr,const vector<vector<double>> & theta_sign,rng<RNGTYPE> & stream)580 void paircompare_theta_update (Matrix<>& theta, const Matrix<>& Ystar,
581 const Matrix<unsigned int>& MD,
582 const Matrix<>& alpha,
583 const Matrix<unsigned int>& theta_n,
584 const Matrix<>& theta_eq,
585 const Matrix<>& theta_ineq,
586 const vector< vector < double* > >& theta_Ystar_ptr,
587 const vector< vector < double* > >& theta_alpha_ptr,
588 const vector< vector < double* > >& theta_theta_ptr,
589 const vector< vector < double > >& theta_sign,
590 rng<RNGTYPE>& stream){
591
592 const unsigned int J = theta.rows();
593 for (unsigned int j = 0; j < J; ++j){
594 double xx = 0.0;
595 double xz = 0.0;
596 for (unsigned int i = 0; i < theta_n[j]; ++i){
597 const double x_ji = theta_sign[j][i] * *theta_alpha_ptr[j][i];
598 const double z_ji = *theta_Ystar_ptr[j][i] + theta_sign[j][i] *
599 *theta_alpha_ptr[j][i] * *theta_theta_ptr[j][i];
600 /*
601 The reason for the above plus sign: -(-theta_sign[j][i])=+theta_sign[j][i].
602 The other theta's sign in a comparison is the opposite to the theta in scrutiny.
603 */
604 xx += x_ji * x_ji;
605 xz += x_ji * z_ji;
606 }
607 // prior for theta is N(0,1), therefore the posterior mean and variance are as follow
608
609 const double v = 1.0/(xx + 1.0);
610 const double m = v * (xz);
611 // no equality constraints
612 if (theta_eq(j) == -999) {
613 if (theta_ineq(j) == 0) { // no inequality constraint
614 theta(j) = stream.rnorm(m, std::sqrt(v));
615 } else if (theta_ineq(j) > 0) { // theta[j] > 0
616 theta(j) = stream.rtbnorm_combo(m, v, 0);
617 } else { // theta[j] < 0
618 theta(j) = stream.rtanorm_combo(m, v, 0);
619 }
620 } else { // equality constraints
621 theta(j) = theta_eq(j);
622 }
623 }
624
625 }
626
627
628
629
630
631
632
633
634
635 // update alpha values in paired comparisons model
636 template <typename RNGTYPE>
paircompare_alpha_update(Matrix<> & alpha,const Matrix<> & Ystar,const Matrix<unsigned int> & MD,const Matrix<> & theta,const double & A0,const double & A0a0,const Matrix<unsigned int> & alpha_n,const vector<vector<double * >> & alpha_Ystar_ptr,const vector<vector<double * >> & alpha_theta1_ptr,const vector<vector<double * >> & alpha_theta2_ptr,rng<RNGTYPE> & stream)637 void paircompare_alpha_update (Matrix<>& alpha,
638 const Matrix<>& Ystar,
639 const Matrix<unsigned int>& MD,
640 const Matrix<>& theta,
641 const double& A0,
642 const double& A0a0,
643 const Matrix<unsigned int>& alpha_n,
644 const vector< vector < double* > >& alpha_Ystar_ptr,
645 const vector< vector < double* > >& alpha_theta1_ptr,
646 const vector< vector < double* > >& alpha_theta2_ptr,
647 rng<RNGTYPE>& stream){
648
649 const unsigned int I = alpha.rows();
650 for (unsigned int i = 0; i < I; ++i){
651 double xx = 0.0;
652 double xz = 0.0;
653 for (unsigned int j = 0; j < alpha_n[i]; ++j){
654 const double x_ji = *alpha_theta1_ptr[i][j] - *alpha_theta2_ptr[i][j];
655 const double z_ji = *alpha_Ystar_ptr[i][j];
656 xx += x_ji * x_ji;
657 xz += x_ji * z_ji;
658 }
659 // prior for alpha is N(0,1), therefore the posterior mean and variance are as follow. Same as Bayesian regression
660 const double v = 1.0/(xx + A0);
661 const double m = v * (xz + A0a0);
662 alpha(i) = stream.rnorm(m, std::sqrt(v));
663 }
664 }
665
666
667
668
669
670
671
672
673
674
675
676
677 // factor analysis model with normal mean 0, precision F0 prior on
678 // factor scores
679 // X follows a multivariate normal distribution
680 // Lambda is the matrix of factor loadings
681 // Psi_inv is the inverse of the uniqueness matrix
682 // N is number of observations
683 // D is the number of factors
684 // this function draws the factor scores
685 //
686 // IMPORTANT
687 // ***********Psi_inv IS ASSUMED TO DIAGONAL ***********
688 template <typename RNGTYPE>
689 void
NormNormfactanal_phi_draw(Matrix<> & phi,const Matrix<> & F0,const Matrix<> & Lambda,const Matrix<> & Psi_inv,const Matrix<> & X,const int & N,const int & D,rng<RNGTYPE> & stream)690 NormNormfactanal_phi_draw(Matrix<> &phi,
691 const Matrix<> &F0,
692 const Matrix<> &Lambda,
693 const Matrix<> &Psi_inv,
694 const Matrix<> &X,
695 const int& N, const int& D,
696 rng <RNGTYPE>& stream){
697 // If Psi_inv is *not* diagonal then use:
698 // Matrix<double> phi_post_var = invpd(F0 + t(Lambda) * Psi_inv *
699 // Lambda);
700 //Instead of the following 2 lines:
701 const Matrix<> AAA = scythe::sqrt(Psi_inv) * Lambda;
702 const Matrix<> phi_post_var = invpd(F0 + crossprod(AAA));
703 const Matrix<> phi_post_C = cholesky(phi_post_var);
704 for (int i=0; i<N; ++i){
705 const Matrix<> phi_post_mean = phi_post_var *
706 (t(Lambda) * Psi_inv * t(X(i,_)));
707 const Matrix<> phi_samp = gaxpy(phi_post_C, stream.rnorm(D, 1, 0.0, 1.0),
708 phi_post_mean);
709 for (int j=0; j<D; ++j)
710 phi(i,j) = phi_samp(j);
711 }
712 }
713
714
715 // samples the Psi matrix for a Normal theory factor model with IG
716 // prior on diag elements of Psi
717 template <typename RNGTYPE>
718 void
NormIGfactanal_Psi_draw(Matrix<> & Psi,const Matrix<> & X,const Matrix<> & phi,const Matrix<> & Lambda,const Matrix<> & a0,const Matrix<> & b0,const int & K,const int & N,rng<RNGTYPE> & stream)719 NormIGfactanal_Psi_draw(Matrix<> &Psi, const Matrix<> &X,
720 const Matrix<> &phi,
721 const Matrix<> &Lambda,
722 const Matrix<> &a0,
723 const Matrix<> &b0,
724 const int& K, const int& N,
725 rng<RNGTYPE>& stream){
726 for (int i=0; i<K; ++i){
727 const Matrix<double> epsilon = gaxpy(phi, -1*(t(Lambda(i,_))), X(_,i));
728 const Matrix<double> SSE = crossprod(epsilon);
729 const double a1 = (a0[i] + N)*0.5;
730 const double b1 = (b0[i] + SSE[0])*0.5;
731 Psi(i,i) = stream.rigamma(a1, b1);
732 }
733 }
734
735
736
737 // Psi_inv assumed diagnonal
738 // this function draws the factor loading matrix
739 template <typename RNGTYPE>
740 void
NormNormfactanal_Lambda_draw(Matrix<> & Lambda,const Matrix<bool> & Lambda_free_indic,const Matrix<> & Lambda_prior_mean,const Matrix<> & Lambda_prior_prec,const Matrix<> & phi,const Matrix<> & X,const Matrix<> & Psi_inv,const Matrix<> & Lambda_ineq,const unsigned int D,const unsigned int K,rng<RNGTYPE> & stream)741 NormNormfactanal_Lambda_draw(Matrix<>& Lambda,
742 const Matrix<bool> &Lambda_free_indic,
743 const Matrix<> &Lambda_prior_mean,
744 const Matrix<> &Lambda_prior_prec,
745 const Matrix<> &phi,
746 const Matrix<> &X,
747 const Matrix<> &Psi_inv,
748 const Matrix<> &Lambda_ineq,
749 const unsigned int D, const unsigned int K,
750 rng<RNGTYPE>& stream)
751 {
752
753 for (unsigned int i = 0; i < K; ++i) {
754 const Matrix<bool> free_indic = t(Lambda_free_indic(i,_));
755 // end replacement
756
757 if (sumc(free_indic)(0) > 0 && sumc(! free_indic)(0) > 0) {
758 // both constrnd & unconstrnd
759 const Matrix<> phifree_i = t(selif(t(phi), free_indic));
760 const Matrix<> mulamfree_i = selif(t(Lambda_prior_mean(i,_)),
761 free_indic); // prior mean
762 const Matrix<> hold = selif(t(Lambda_prior_prec(i,_)),
763 free_indic);
764 Matrix<> sig2lamfree_inv_i
765 = eye<double>(hold.rows()); // prior prec
766
767 for (unsigned int j = 0; j < (hold.rows()); ++j)
768 sig2lamfree_inv_i(j,j) = hold[j];
769
770 const Matrix<> Lambdacon_i = selif(t(Lambda(i,_)), ! free_indic);
771 const Matrix<> phicon_i = t(selif(t(phi), ! free_indic));
772 const Matrix<> newX_i = gaxpy((-1.0*phicon_i), Lambdacon_i,
773 X(_,i));
774 const Matrix<> Lam_post_var = invpd(sig2lamfree_inv_i +
775 Psi_inv(i,i) * crossprod(phifree_i));
776 const Matrix<> Lam_post_C = cholesky(Lam_post_var);
777 const Matrix<> Lam_post_mean = Lam_post_var *
778 (sig2lamfree_inv_i * mulamfree_i + Psi_inv(i,i) *
779 t(phifree_i) * newX_i);
780
781
782 Matrix<> Lambdafree_i =
783 gaxpy(Lam_post_C, stream.rnorm(hold.rows(), 1, 0, 1),Lam_post_mean);
784
785
786
787 // check to see if inequality constraints hold
788 const Matrix<> Lambda_ineq_vec = Lambda_ineq(i,_);
789
790 double ineq_holds = 0;
791 int Lam_count = 0;
792 for (unsigned int j = 0; j < D; ++j) {
793 if (free_indic(j)){
794 ineq_holds = std::min(ineq_holds, Lambda_ineq_vec(j) *
795 Lambdafree_i(Lam_count));
796 ++Lam_count;
797 }
798 }
799
800 while (ineq_holds < 0) {
801 Lambdafree_i = gaxpy(Lam_post_C, stream.rnorm(hold.rows(), 1, 0, 1),
802 Lam_post_mean);
803 Lam_count = 0;
804 double test = 0;
805 for (unsigned int j = 0; j < D; ++j) {
806 if (free_indic(j) == 1) {
807 Matrix<> prodcheck = Lambda_ineq_vec(j)
808 * Lambdafree_i(Lam_count);
809 test = std::min(test, prodcheck(0));
810 ++Lam_count;
811 }
812 }
813 ineq_holds = test;
814 }
815
816 // put draw into Lambda
817 Lam_count = 0;
818 for (unsigned int j = 0; j < D; ++j) {
819 if (free_indic(j) == 1) {
820 Lambda(i,j) = Lambdafree_i(Lam_count);
821 ++Lam_count;
822 }
823 }
824 } else if (sumc(free_indic)(0) > 0) { // just unconstrained
825 const Matrix<> phifree_i = t(selif(t(phi), free_indic));
826 const Matrix<> mulamfree_i = selif(t(Lambda_prior_mean(i,_)),
827 free_indic); // prior mean
828 const Matrix<> hold = selif(t(Lambda_prior_prec(i,_)), free_indic);
829 Matrix<> sig2lamfree_inv_i = eye<double>(hold.rows()); // prior prec
830
831 for (unsigned int j = 0; j < hold.rows(); ++j)
832 sig2lamfree_inv_i(j,j) = hold(j);
833
834 const Matrix<> Lam_post_var = invpd(sig2lamfree_inv_i +
835 Psi_inv(i,i) * crossprod(phifree_i));
836 const Matrix<> Lam_post_C = cholesky(Lam_post_var);
837 const Matrix<> Lam_post_mean = Lam_post_var * (sig2lamfree_inv_i
838 * mulamfree_i +
839 Psi_inv(i,i) *
840 t(phifree_i) * X(_,i));
841 Matrix<> Lambdafree_i = gaxpy(Lam_post_C,
842 stream.rnorm(hold.rows(), 1, 0, 1),
843 Lam_post_mean);
844
845
846 // check to see if inequality constraints hold
847 Matrix<> Lambda_ineq_vec = Lambda_ineq(i,_);
848 double ineq_holds = 0;
849 for (unsigned int j = 0; j < D; ++j) {
850 ineq_holds = std::min(ineq_holds, Lambda_ineq_vec(j)
851 * Lambdafree_i(j));
852 }
853
854 while (ineq_holds < 0) {
855 Lambdafree_i = gaxpy(Lam_post_C, stream.rnorm(hold.rows(), 1, 0, 1),
856 Lam_post_mean);
857 double test = 0;
858 for (unsigned int j = 0; j < D; ++j) {
859 //if (free_indic[j]==1)
860 double prodcheck = Lambda_ineq_vec[j]*Lambdafree_i[j];
861 test = std::min(test, prodcheck);
862 }
863
864 ineq_holds = test;
865 }
866
867 // put draw into Lambda
868 for (unsigned int j = 0; j < D; ++j) {
869 Lambda(i,j) = Lambdafree_i(j);
870 }
871 }
872 }
873 // return(Lambda);
874 }
875
876
877 // update ability parameters (ideal points) for one dimensional
878 // Hierarchical item response model.
879 // 2008-11-18 now deals with PX alpha and holds on to thetahat
880 template <typename RNGTYPE>
hirt_theta_update1(Matrix<> & theta,Matrix<> & thetahat,const Matrix<> & Z,const Matrix<> & eta,const Matrix<> & beta,const Matrix<> & Xj,const double & sigma2,const double & alpha,rng<RNGTYPE> & stream)881 void hirt_theta_update1 (Matrix<>& theta, Matrix<>& thetahat,
882 const Matrix<>& Z,
883 const Matrix<>& eta,
884 const Matrix<>& beta, const Matrix<>& Xj,
885 const double& sigma2,
886 const double& alpha,
887 rng<RNGTYPE>& stream)
888 {
889
890 const unsigned int J = Z.rows();
891 const unsigned int K = Z.cols();
892 // Get level1 prior mean
893 const Matrix<double> Xbeta = (Xj * beta);
894
895 // a and b are backwards here, different parameterizations
896 // common for edu people and us.
897 const Matrix<> b = eta(_, 0); // location
898 const Matrix<> a = eta(_, 1); // relevance
899
900 // calculate the posterior variance outside the justice specific loop
901 const double sig2_inv = 1.0 / sigma2;
902 const Matrix<double> apa = crossprod(a);
903 const Matrix<double> theta_post_var = scythe::invpd(apa + sig2_inv);
904 const double theta_post_sd = std::sqrt(theta_post_var[0]);
905 // sample for each justice
906 for (unsigned int j = 0; j < J; ++j) {
907 thetahat(j) = 0.0;
908 for (unsigned int k = 0; k < K; ++k) {
909 thetahat(j) += a[k] * ( Z(j,k) + b[k]); // bill contribution
910 }
911 thetahat(j) += (Xbeta[j] / sigma2); // j prior level1 contribution
912
913 thetahat(j) *= theta_post_var[0];
914 const double t = thetahat(j) / alpha;
915 theta(j) = stream.rnorm( t , theta_post_sd );
916 }
917 }
918
919
920 // update item (case, roll call) parameters for item response model
921 // note: works only for one-dimensional case
922 // updated for PX (alpha) and hold on to etahat 2008-11-18
923 template <typename RNGTYPE>
hirt_eta_update1(Matrix<> & eta,Matrix<> & etahat,const Matrix<> & Z,const Matrix<> & theta,const Matrix<> & AB0,const Matrix<> & AB0ab0,const double & alpha,rng<RNGTYPE> & stream)924 void hirt_eta_update1 (Matrix<>& eta, Matrix<>& etahat, const Matrix<>& Z,
925 const Matrix<>& theta, const Matrix<>& AB0,
926 const Matrix<>& AB0ab0,
927 const double& alpha,
928 rng<RNGTYPE>& stream)
929 {
930
931 // define constants
932 const unsigned int J = theta.rows();
933 const unsigned int K = Z.cols();
934
935 // perform update
936 //const Matrix<double> Ttheta_star = t(cbind(-1.0*ones<double>(J,1),theta)); // only needed for option 2
937 const Matrix<> tpt(2,2);
938 for (unsigned int i = 0; i < J; ++i) {
939 const double theta_i = theta(i);
940 tpt(0,1) -= theta_i;
941 tpt(1,1) += std::pow(theta_i, 2.0);
942 }
943 tpt(1,0) = tpt(0,1);
944 tpt(0,0) = J;
945 const Matrix<> eta_post_var = invpd(tpt + AB0);
946 const Matrix<> eta_post_C = cholesky(eta_post_var);
947
948 for (unsigned int k = 0; k < K; ++k) {
949 const Matrix<> TZ(2, 1);
950 for (unsigned int j = 0; j < J; ++j) {
951 TZ[0] -= Z(j,k);
952 TZ[1] += Z(j,k) * theta[j];
953 }
954 Matrix<> eta_post_mean = eta_post_var * (TZ + AB0ab0);
955 etahat(k,0) = eta_post_mean(0);
956 etahat(k,1) = eta_post_mean(1);
957 eta_post_mean /= alpha;
958 const Matrix<> new_eta = gaxpy(eta_post_C, stream.rnorm(2, 1, 0, 1),
959 (eta_post_mean) );
960 eta(k,0) = new_eta(0);
961 eta(k,1) = new_eta(1);
962
963 //Rprintf("\n\a: %3.1f,b:%3.1f ",eta(k,0),eta(k,1));
964
965 }
966 }
967
968
969
970
971
972
973
974 #endif
975
976
977
978