1 //////////////////////////////////////////////////////////////////// 2 // cMCMChregress.cc 3 // 4 // cMCMChregress samples from the posterior distribution of a 5 // Gaussian hierarchical linear regression model 6 // 7 // The code uses Algorithm 2 of Chib & Carlin (1999) for efficient 8 // inference of (\beta | Y, sigma^2, Vb). 9 // 10 // Chib, S. & Carlin, B. P. (1999) On MCMC sampling in hierarchical 11 // longitudinal models, Statistics and Computing, 9, 17-26 12 // 13 //////////////////////////////////////////////////////////////////// 14 // 15 // Original code by Ghislain Vieilledent, may 2011 16 // CIRAD UR B&SEF 17 // ghislain.vieilledent@cirad.fr / ghislainv@gmail.com 18 // 19 //////////////////////////////////////////////////////////////////// 20 // 21 // The initial version of this file was generated by the 22 // auto.Scythe.call() function in the MCMCpack R package 23 // written by: 24 // 25 // Andrew D. Martin 26 // Dept. of Political Science 27 // Washington University in St. Louis 28 // admartin@wustl.edu 29 // 30 // Kevin M. Quinn 31 // Dept. of Government 32 // Harvard University 33 // kevin_quinn@harvard.edu 34 // 35 // This software is distributed under the terms of the GNU GENERAL 36 // PUBLIC LICENSE Version 2, June 1991. See the package LICENSE 37 // file for more information. 38 // 39 // Copyright (C) 2011 Andrew D. Martin and Kevin M. Quinn 40 // 41 //////////////////////////////////////////////////////////////////// 42 // 43 // Revisions: 44 // - This file was initially generated on Wed May 4 10:42:50 2011 45 // - G. Vieilledent, on May 4 2011 46 // 47 //////////////////////////////////////////////////////////////////// 48 49 #include "matrix.h" 50 #include "distributions.h" 51 #include "stat.h" 52 #include "la.h" 53 #include "ide.h" 54 #include "smath.h" 55 #include "MCMCrng.h" 56 #include "MCMCfcds.h" 57 58 #include <R.h> // needed to use Rprintf() 59 #include <R_ext/Utils.h> // needed to allow user interrupts 60 61 using namespace scythe; 62 using namespace std; 63 64 extern "C"{ 65 66 /* Gibbs sampler function */ cMCMChregress(const int * ngibbs,const int * nthin,const int * nburn,const int * nobs,const int * ngroup,const int * np,const int * nq,const int * IdentGroup,const double * Y_vect,const double * X_vect,const double * W_vect,double * beta_vect,double * b_vect,double * Vb_vect,double * V,const double * mubeta_vect,const double * Vbeta_vect,const double * r,const double * R_vect,const double * s1_V,const double * s2_V,double * Deviance,double * Y_pred,const int * seed,const int * verbose)67 void cMCMChregress ( 68 69 // Constants and data 70 const int *ngibbs, const int *nthin, const int *nburn, // Number of iterations, burning and samples 71 const int *nobs, const int *ngroup, // Constants 72 const int *np, const int *nq, // Number of fixed and random covariates 73 const int *IdentGroup, // Vector of group 74 const double *Y_vect, // Observed response variable 75 const double *X_vect, // Covariate for fixed effects 76 const double *W_vect, // Covariate for random effects 77 // Parameters to save 78 double *beta_vect, // Fixed effects 79 double *b_vect, // Random effects 80 double *Vb_vect, // Variance of random effects 81 double *V, // Variance of residuals 82 // Defining priors 83 const double *mubeta_vect, const double *Vbeta_vect, 84 const double *r, const double *R_vect, 85 const double *s1_V, const double *s2_V, 86 // Diagnostic 87 double *Deviance, 88 double *Y_pred, // Fitted values (predictive posterior mean) 89 // Seeds 90 const int *seed, 91 // Verbose 92 const int *verbose 93 94 ) { 95 96 //////////////////////////////////////////////////////////////////////////////// 97 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 98 // Defining and initializing objects 99 100 /////////////////////////// 101 // Redefining constants // 102 const int NGIBBS=ngibbs[0]; 103 const int NTHIN=nthin[0]; 104 const int NBURN=nburn[0]; 105 const int NSAMP=(NGIBBS-NBURN)/NTHIN; 106 const int NOBS=nobs[0]; 107 const int NGROUP=ngroup[0]; 108 const int NP=np[0]; 109 const int NQ=nq[0]; 110 111 /////////////// 112 // Constants // 113 // Number of observations by group k 114 int *nobsk = new int[NGROUP]; 115 for (int k=0; k<NGROUP; k++) { 116 nobsk[k]=0; 117 for (int n=0; n<NOBS; n++) { 118 if (IdentGroup[n]==k) { 119 nobsk[k]+=1; 120 } 121 } 122 } 123 // Position of each group in the data-set 124 int **posk_arr = new int*[NGROUP]; 125 for (int k=0; k<NGROUP; k++) { 126 posk_arr[k] = new int[nobsk[k]]; 127 int repk=0; 128 for (int n=0; n<NOBS; n++) { 129 if (IdentGroup[n]==k) { 130 posk_arr[k][repk]=n; 131 repk++; 132 } 133 } 134 } 135 // Small fixed matrices indexed on k for data access 136 Matrix<double> *Yk_arr = new Matrix<double>[NGROUP]; 137 Matrix<double> *Xk_arr = new Matrix<double>[NGROUP]; 138 Matrix<double> *Wk_arr = new Matrix<double>[NGROUP]; 139 for(int k=0; k<NGROUP; k++) { 140 Xk_arr[k] = Matrix<double>(nobsk[k],NP); 141 Wk_arr[k] = Matrix<double>(nobsk[k],NQ); 142 Yk_arr[k] = Matrix<double>(nobsk[k],1); 143 for (int m=0; m<nobsk[k]; m++) { 144 for (int p=0; p<NP; p++) { 145 Xk_arr[k](m,p)=X_vect[p*NOBS+posk_arr[k][m]]; 146 } 147 for (int q=0; q<NQ; q++) { 148 Wk_arr[k](m,q)=W_vect[q*NOBS+posk_arr[k][m]]; 149 } 150 Yk_arr[k](m,0)=Y_vect[posk_arr[k][m]]; 151 Y_pred[posk_arr[k][m]]=0; // We initialize Y_pred to zero to compute the predictive posterior mean 152 } 153 } 154 Matrix<double> *tXk_arr = new Matrix<double>[NGROUP]; 155 Matrix<double> *tWk_arr = new Matrix<double>[NGROUP]; 156 Matrix<double> *cpXk_arr = new Matrix<double>[NGROUP]; 157 Matrix<double> *tXWk_arr = new Matrix<double>[NGROUP]; 158 Matrix<double> *tWXk_arr = new Matrix<double>[NGROUP]; 159 Matrix<double> *tXYk_arr = new Matrix<double>[NGROUP]; 160 Matrix<double> *tWYk_arr = new Matrix<double>[NGROUP]; 161 for(int k=0; k<NGROUP; k++) { 162 tXk_arr[k] = t(Xk_arr[k]); 163 tWk_arr[k] = t(Wk_arr[k]); 164 cpXk_arr[k] = crossprod(Xk_arr[k]); 165 tXWk_arr[k] = t(Xk_arr[k])*Wk_arr[k]; 166 tWXk_arr[k] = t(Wk_arr[k])*Xk_arr[k]; 167 tXYk_arr[k] = t(Xk_arr[k])*Yk_arr[k]; 168 tWYk_arr[k] = t(Wk_arr[k])*Yk_arr[k]; 169 } 170 171 //////////////////// 172 // Priors objects // 173 Matrix<double> mubeta(NP,1,mubeta_vect); 174 Matrix<double> Vbeta(NP,NP,Vbeta_vect); 175 Matrix<double> R(NQ,NQ,R_vect); 176 177 ///////////////////////////////////// 178 // Initializing running parameters // 179 Matrix<double> *bk_run = new Matrix<double>[NGROUP]; // Random effects 180 for (int k=0;k<NGROUP;k++) { 181 bk_run[k] = Matrix<double>(NQ,1); 182 for (int q=0; q<NQ; q++) { 183 bk_run[k](q)=b_vect[q*NGROUP*NSAMP+k*NSAMP]; 184 } 185 } 186 Matrix<double> beta_run(NP,1,false); // Unicolumn matrix of fixed effects 187 for (int p=0; p<NP; p++) { 188 beta_run(p)=beta_vect[p*NSAMP]; 189 } 190 Matrix<double> Vb_run(NQ,NQ,true,0.0); 191 for (int q=0; q<NQ; q++) { 192 for (int qprim=0; qprim<NQ; qprim++) { 193 Vb_run(q,qprim)=Vb_vect[qprim*NQ*NSAMP+q*NSAMP]; 194 } 195 } 196 double V_run=V[0]; 197 double Deviance_run=Deviance[0]; 198 199 //////////// 200 // Message// 201 Rprintf("\nRunning the Gibbs sampler. It may be long, keep cool :)\n\n"); 202 R_FlushConsole(); 203 //R_ProcessEvents(); for windows 204 205 ///////////////////// 206 // Set random seed // 207 mersenne myrng; 208 myrng.initialize(*seed); 209 210 211 /////////////////////////////////////////////////////////////////////////////////////// 212 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 213 // Gibbs sampler (see Chib et Carlin, 1999 p.4 "Blocking for Gaussian mixed models") 214 215 for (int g=0;g<NGIBBS;g++) { 216 217 ////////////////////////////////// 218 // vector beta: Gibbs algorithm // 219 220 // invVi, sum_V and sum_v 221 // see http://en.wikipedia.org/wiki/Woodbury_matrix_identity with A=(1/V_run)*Id 222 Matrix<double> sum_V(NP,NP); 223 Matrix<double> sum_v(NP,1); 224 for (int k=0; k<NGROUP; k++) { 225 sum_V += (1/V_run)*cpXk_arr[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWXk_arr[k]; 226 sum_v += (1/V_run)*tXYk_arr[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWYk_arr[k]; 227 } 228 229 // big_V 230 Matrix<double> big_V=invpd(invpd(Vbeta)+sum_V/V_run); 231 232 // small_v 233 Matrix<double> small_v=invpd(Vbeta)*mubeta+sum_v/V_run; 234 235 // Draw in the posterior distribution 236 beta_run=myrng.rmvnorm(big_V*small_v,big_V); 237 238 239 /////////////////////////////// 240 // vector b: Gibbs algorithm // 241 242 // Loop on group 243 for (int k=0; k<NGROUP; k++) { 244 245 // big_Vk 246 Matrix<double> big_Vk=invpd(invpd(Vb_run)+crossprod(Wk_arr[k])/V_run); 247 248 // small_vk 249 Matrix<double> small_vk=(t(Wk_arr[k])*(Yk_arr[k]-Xk_arr[k]*beta_run))/V_run; 250 251 // Draw in the posterior distribution 252 bk_run[k]=myrng.rmvnorm(big_Vk*small_vk,big_Vk); 253 } 254 255 256 //////////////////////////////////////////// 257 // vector of variance Vb: Gibbs algorithm // 258 Matrix<double> SSBb(NQ,NQ,true,0.0); 259 for(int k=0; k<NGROUP; k++) { 260 SSBb+=bk_run[k]*t(bk_run[k]); 261 } 262 int Vb_dof=(*r)+(NGROUP); 263 Matrix<double> Vb_scale = invpd(SSBb+(*r)*R); 264 Vb_run=invpd(myrng.rwish(Vb_dof,Vb_scale)); 265 266 267 //////////////// 268 // variance V // 269 270 // e 271 Matrix<double> e(1,1,true,0.0); 272 for (int k=0; k<NGROUP; k++) { 273 e+=crossprod(Yk_arr[k]-Xk_arr[k]*beta_run-Wk_arr[k]*bk_run[k]); 274 } 275 276 // Parameters 277 double S1=*s1_V+(NOBS/2); //shape 278 double S2=*s2_V+0.5*e(0); //rate 279 280 // Draw in the posterior distribution 281 V_run=1/myrng.rgamma(S1,S2); 282 283 ////////////////////////////////////////////////// 284 //// Deviance 285 286 // logLikelihood 287 double logLk=0; 288 for (int k=0; k<NGROUP; k++) { 289 for (int m=0; m<nobsk[k]; m++) { 290 // Which one ? 291 int w=posk_arr[k][m]; 292 // Y_hat 293 Matrix<double> Y_hat=Xk_arr[k](m,_)*beta_run+Wk_arr[k](m,_)*bk_run[k]; 294 // L 295 logLk+=log(dnorm(Y_vect[w],Y_hat(0),sqrt(V_run))); 296 } 297 } 298 299 // Deviance 300 Deviance_run=-2*logLk; 301 302 303 ////////////////////////////////////////////////// 304 // Output 305 if(((g+1)>NBURN) && (((g+1)%(NTHIN))==0)) { 306 int isamp=((g+1)-NBURN)/(NTHIN); 307 for (int p=0; p<NP; p++) { 308 beta_vect[p*NSAMP+(isamp-1)]=beta_run(p); 309 } 310 for (int k=0; k<NGROUP; k++) { 311 for (int q=0; q<NQ; q++) { 312 b_vect[q*NGROUP*NSAMP+k*NSAMP+(isamp-1)]=bk_run[k](q); 313 } 314 } 315 for (int q=0; q<NQ; q++) { 316 for (int qprim=0; qprim<NQ; qprim++) { 317 Vb_vect[qprim*NQ*NSAMP+q*NSAMP+(isamp-1)]=Vb_run(q,qprim); 318 } 319 } 320 V[isamp-1]=V_run; 321 Deviance[isamp-1]=Deviance_run; 322 for (int k=0; k<NGROUP; k++) { 323 for (int m=0; m<nobsk[k]; m++) { 324 // Which one ? 325 int w=posk_arr[k][m]; 326 // Y_hat 327 Matrix<double> Y_hat=Xk_arr[k](m,_)*beta_run+Wk_arr[k](m,_)*bk_run[k]; 328 // Y_pred 329 Y_pred[w]+=Y_hat(0)/NSAMP; 330 } 331 } 332 } 333 334 335 ////////////////////////////////////////////////// 336 // Progress bar 337 double Perc=100*(g+1)/(NGIBBS); 338 if(((g+1)%(NGIBBS/100))==0 && (*verbose==1)){ 339 Rprintf("*"); 340 R_FlushConsole(); 341 //R_ProcessEvents(); for windows 342 if(((g+1)%(NGIBBS/10))==0){ 343 Rprintf(":%.1f%%\n",Perc); 344 R_FlushConsole(); 345 //R_ProcessEvents(); for windows 346 } 347 } 348 349 350 ////////////////////////////////////////////////// 351 // User interrupt 352 R_CheckUserInterrupt(); // allow user interrupts 353 354 } // end MCMC loop 355 356 /////////////// 357 // Delete memory allocation 358 delete[] nobsk; 359 for(int k=0; k<NGROUP; k++) { 360 delete[] posk_arr[k]; 361 } 362 delete[] posk_arr; 363 delete[] Yk_arr; 364 delete[] Xk_arr; 365 delete[] Wk_arr; 366 delete[] tXk_arr; 367 delete[] tWk_arr; 368 delete[] cpXk_arr; 369 delete[] tXWk_arr; 370 delete[] tWXk_arr; 371 delete[] tXYk_arr; 372 delete[] tWYk_arr; 373 delete[] bk_run; 374 375 } // end cMCMChregress 376 377 } // end extern "C" 378 379 //////////////////////////////////////////////////////////////////// 380 // END 381 //////////////////////////////////////////////////////////////////// 382