1 //////////////////////////////////////////////////////////////////////////
2 // cMCMCfactanal.cc is C++ code to estimate a factor analysis model
3 //
4 // Andrew D. Martin
5 // Dept. of Political Science
6 // Washington University in St. Louis
7 // admartin@wustl.edu
8 //
9 // Kevin M. Quinn
10 // Dept. of Government
11 // Harvard University
12 // kevin_quinn@harvard.edu
13 //
14 // This software is distributed under the terms of the GNU GENERAL
15 // PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
16 // file for more information.
17 //
18 // revised version of older MCMCfactanal 5/11/2004 KQ
19 // updated to new verion of scythe 7/25/2004 ADM
20 // updated to Scythe 1.0.X 7/10/2007 ADM
21 // finished update to Scythe 1.0.X 7/30/2007 KQ
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 #ifndef CMCMCFACTANAL_CC
30 #define CMCMCFACTANAL_CC
31
32 #include "matrix.h"
33 #include "algorithm.h"
34 #include "distributions.h"
35 #include "stat.h"
36 #include "la.h"
37 #include "ide.h"
38 #include "smath.h"
39 #include "MCMCrng.h"
40 #include "MCMCfcds.h"
41
42 #include <R.h> // needed to use Rprintf()
43 #include <R_ext/Utils.h> // needed to allow user interrupts
44
45 typedef Matrix<double,Row,View> rmview;
46
47 using namespace std;
48 using namespace scythe;
49
50 template <typename RNGTYPE>
MCMCfactanal_impl(rng<RNGTYPE> & stream,const Matrix<> & X,Matrix<> & Lambda,Matrix<> & Psi,Matrix<> & Psi_inv,const Matrix<> & Lambda_eq,const Matrix<> & Lambda_ineq,const Matrix<> & Lambda_prior_mean,const Matrix<> & Lambda_prior_prec,const Matrix<> & a0,const Matrix<> & b0,unsigned int burnin,unsigned int mcmc,unsigned int thin,unsigned int verbose,unsigned int storescores,Matrix<> & result)51 void MCMCfactanal_impl (rng<RNGTYPE>& stream, const Matrix<>& X,
52 Matrix<>& Lambda,
53 Matrix<>& Psi, Matrix<>& Psi_inv,
54 const Matrix<>& Lambda_eq,
55 const Matrix<>& Lambda_ineq,
56 const Matrix<>& Lambda_prior_mean,
57 const Matrix<>& Lambda_prior_prec,
58 const Matrix<>& a0, const Matrix<>& b0,
59 unsigned int burnin, unsigned int mcmc,
60 unsigned int thin,
61 unsigned int verbose,
62 unsigned int storescores, Matrix<>& result) {
63
64 // constants
65 const unsigned int K = X.cols(); // number of manifest variables
66 const unsigned int N = X.rows(); // number of observations
67 const unsigned int D = Lambda.cols(); // number of factors
68 const unsigned int tot_iter = burnin + mcmc;
69 const unsigned int nsamp = mcmc / thin;
70 const Matrix<> I = eye(D);
71 const Matrix<> Lambda_free_indic = Matrix<>(K, D);
72 for (unsigned int i=0; i<(K*D); ++i){
73 if (Lambda_eq[i] == -999) Lambda_free_indic[i] = 1.0;
74 }
75
76
77 // starting value for phi
78 Matrix<> phi = Matrix<>(N,D);
79
80 // storage matrices (row major order)
81 Matrix<> Lambda_store = Matrix<>(nsamp, K*D);
82 Matrix<> Psi_store = Matrix<>(nsamp, K);
83 Matrix<> phi_store;
84 if (storescores==1){
85 phi_store = Matrix<double>(nsamp, N*D);
86 }
87
88
89 unsigned int count = 0;
90 // sampling begins here
91 for (unsigned int iter=0; iter < tot_iter; ++iter){
92
93 // sample phi
94 NormNormfactanal_phi_draw(phi, I, Lambda, Psi_inv,
95 X, N, D, stream);
96
97
98 // sample Lambda
99 NormNormfactanal_Lambda_draw(Lambda, Lambda_free_indic,
100 Lambda_prior_mean,
101 Lambda_prior_prec,
102 phi, X, Psi_inv,
103 Lambda_ineq,
104 D, K, stream);
105
106
107 // sample Psi
108 NormIGfactanal_Psi_draw(Psi, X, phi, Lambda,
109 a0, b0, K, N, stream);
110
111
112 for (unsigned int i=0; i<K; ++i)
113 Psi_inv(i,i) = 1.0 / Psi(i,i);
114
115
116
117 // print results to screen
118 if (verbose > 0 && iter % verbose == 0){
119 Rprintf("\n\nMCMCfactanal iteration %i of %i \n", (iter+1), tot_iter);
120 Rprintf("Lambda = \n");
121 for (unsigned int i=0; i<K; ++i){
122 for (unsigned int j=0; j<D; ++j){
123 Rprintf("%10.5f", Lambda(i,j));
124 }
125 Rprintf("\n");
126 }
127 Rprintf("diag(Psi) = \n");
128 for (unsigned int i=0; i<K; ++i){
129 Rprintf("%10.5f", Psi(i,i));
130 }
131 Rprintf("\n");
132 }
133
134 // store results
135 if ((iter % thin)==0 && iter >= burnin ) {
136 // store Lambda
137 //Matrix<> Lambda_store_vec = Lambda.resize(1, K*D, true);
138 //for (int l=0; l<K*D; ++l)
139 //Lambda_store(count, l) = Lambda(l);
140 rmview(Lambda_store(count, _)) = Lambda;
141
142 // store Psi
143 for (unsigned int i=0; i<K; ++i)
144 Psi_store(count, i) = Psi(i,i);
145 // stop phi
146 if (storescores==1){
147 //Matrix<> phi_store_vec = phi.resize(1, N*D, true);
148 //for (int l=0; l<N*D; ++l)
149 // phi_store(count, l) = phi(l);
150 rmview(phi_store(count, _)) = phi;
151 }
152 count++;
153 }
154
155 // allow user interrupts
156 R_CheckUserInterrupt();
157
158 } // end Gibbs loop
159
160 // return output
161 //Matrix<double> result;
162 result = cbind(Lambda_store, Psi_store);
163 if(storescores == 1) {
164 result = cbind(result, phi_store);
165 }
166
167 }
168
169
170 extern "C" {
171
cMCMCfactanal(double * sampledata,const int * samplerow,const int * samplecol,const double * Xdata,const int * Xrow,const int * Xcol,const int * burnin,const int * mcmc,const int * thin,const int * uselecuyer,const int * seedarray,const int * lecuyerstream,const int * verbose,const double * Lambdadata,const int * Lambdarow,const int * Lambdacol,const double * Psidata,const int * Psirow,const int * Psicol,const double * Lameqdata,const int * Lameqrow,const int * Lameqcol,const double * Lamineqdata,const int * Lamineqrow,const int * Lamineqcol,const double * Lampmeandata,const int * Lampmeanrow,const int * Lampmeancol,const double * Lampprecdata,const int * Lampprecrow,const int * Lamppreccol,const double * a0data,const int * a0row,const int * a0col,const double * b0data,const int * b0row,const int * b0col,const int * storescores)172 void cMCMCfactanal(double *sampledata, const int *samplerow,
173 const int *samplecol, const double *Xdata,
174 const int *Xrow, const int *Xcol, const int *burnin,
175 const int *mcmc, const int *thin, const int *uselecuyer,
176 const int *seedarray, const int *lecuyerstream,
177 const int *verbose, const double *Lambdadata,
178 const int *Lambdarow, const int *Lambdacol,
179 const double *Psidata, const int *Psirow,
180 const int *Psicol, const double *Lameqdata,
181 const int *Lameqrow, const int *Lameqcol,
182 const double *Lamineqdata, const int *Lamineqrow,
183 const int *Lamineqcol, const double *Lampmeandata,
184 const int *Lampmeanrow, const int *Lampmeancol,
185 const double *Lampprecdata, const int *Lampprecrow,
186 const int *Lamppreccol, const double *a0data,
187 const int *a0row, const int *a0col,
188 const double *b0data, const int *b0row,
189 const int *b0col, const int *storescores) {
190
191 // pull together Matrix objects
192 const Matrix <> X(*Xrow, *Xcol, Xdata);
193 Matrix <> Lambda(*Lambdarow, *Lambdacol, Lambdadata);
194 Matrix <> Psi(*Psirow, *Psicol, Psidata);
195 Matrix <> Psi_inv = invpd(Psi);
196 const Matrix <> Lambda_eq(*Lameqrow, *Lameqcol,
197 Lameqdata);
198 const Matrix <> Lambda_ineq(*Lamineqrow, *Lamineqcol,
199 Lamineqdata);
200 const Matrix <> Lambda_prior_mean(*Lampmeanrow,
201 *Lampmeancol,
202 Lampmeandata);
203 const Matrix <> Lambda_prior_prec(*Lampprecrow,
204 *Lamppreccol,
205 Lampprecdata);
206 const Matrix <> a0(*a0row, *a0col, a0data);
207 const Matrix <> b0(*b0row, *b0col, b0data);
208
209 Matrix<> storagematrix;
210
211 MCMCPACK_PASSRNG2MODEL(MCMCfactanal_impl, X, Lambda, Psi, Psi_inv,
212 Lambda_eq, Lambda_ineq, Lambda_prior_mean,
213 Lambda_prior_prec,
214 a0, b0, *burnin, *mcmc, *thin,
215 *verbose, *storescores, storagematrix);
216
217 const unsigned int size = *samplerow * *samplecol;
218 for (unsigned int i = 0; i < size; ++i)
219 sampledata[i] = storagematrix(i);
220 }
221 }
222
223 #endif
224
225