1 //////////////////////////////////////////////////////////////////////////
2 // MCMCpaircompare.cc is C++ code to estimate a pairwise comparison model.
3 //
4 // KQ 3/18/2015
5 //
6 // Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
7 // Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
8 // and Jong Hee Park
9 //////////////////////////////////////////////////////////////////////////
10
11 #ifndef MCMCPAIRCOMPARE_CC
12 #define MCMCPAIRCOMPARE_CC
13
14 #include<vector>
15
16 #include "MCMCrng.h"
17 #include "MCMCfcds.h"
18 #include "matrix.h"
19 #include "distributions.h"
20 #include "stat.h"
21 #include "la.h"
22 #include "ide.h"
23 #include "smath.h"
24 #include "rng.h"
25 #include "mersenne.h"
26 #include "lecuyer.h"
27
28 #include <R.h> // needed to use Rprintf()
29 #include <R_ext/Utils.h> // needed to allow user interrupts
30
31 using namespace std;
32 using namespace scythe;
33
34 /* MCMCpaircompare implementation. */
35 template <typename RNGTYPE>
MCMCpaircompare_impl(rng<RNGTYPE> & stream,const Matrix<unsigned int> & MD,Matrix<> & theta,Matrix<> & alpha,const Matrix<> & theta_eq,const Matrix<> & theta_ineq,const double a0,const double A0,const int alphafixed,const unsigned int burnin,const unsigned int mcmc,const unsigned int thin,const unsigned int verbose,const bool storealpha,const bool storetheta,double * sampledata,const unsigned int samplesize)36 void MCMCpaircompare_impl (rng<RNGTYPE>& stream,
37 const Matrix<unsigned int>& MD,
38 Matrix<>& theta,
39 Matrix<>& alpha,
40 const Matrix<>& theta_eq,
41 const Matrix<>& theta_ineq,
42 const double a0,
43 const double A0,
44 const int alphafixed,
45 const unsigned int burnin,
46 const unsigned int mcmc,
47 const unsigned int thin,
48 const unsigned int verbose,
49 const bool storealpha,
50 const bool storetheta,
51 double* sampledata,
52 const unsigned int samplesize){
53
54
55 // constants
56 const unsigned int N = MD.rows();
57 const unsigned int J = theta.rows();
58 const unsigned int I = alpha.rows();
59 const unsigned int tot_iter = burnin + mcmc;
60 const unsigned int nsamp = mcmc / thin;
61
62
63 // starting values for Ystar
64 Matrix<> Ystar(N,1);
65
66
67 // pre-compute what can be pre-computed
68 const double A0a0 = A0 * a0;
69 Matrix<unsigned int> theta_n(J,1,true,0); // vector of 0s. Item j's total instances of being compared.
70 Matrix<unsigned int> alpha_n(I,1,true,0); // vector of 0s. Judge i's total instances of making comparisons.
71 for (unsigned int i = 0; i < N; ++i){
72 alpha_n(MD(i,0)) += 1;
73 theta_n(MD(i,1)) += 1;
74 theta_n(MD(i,2)) += 1;
75 }
76
77 vector< vector< double* > > theta_Ystar_ptr;
78 vector< vector< double* > > theta_alpha_ptr;
79 vector< vector< double* > > theta_theta_ptr;
80 vector< vector< double > > theta_sign;
81
82 vector< vector< double* > > alpha_Ystar_ptr;
83 vector< vector< double* > > alpha_theta1_ptr;
84 vector< vector< double* > > alpha_theta2_ptr;
85
86 theta_Ystar_ptr.reserve(J);
87 theta_alpha_ptr.reserve(J);
88 theta_theta_ptr.reserve(J);
89 theta_sign.reserve(J);
90
91 alpha_Ystar_ptr.reserve(I);
92 alpha_theta1_ptr.reserve(I);
93 alpha_theta2_ptr.reserve(I);
94
95 for (unsigned int j = 0; j < J; ++j){
96 vector< double* > Ystar_j_ptr;
97 vector< double* > alpha_j_ptr;
98 vector< double* > theta_j_ptr;
99 vector< double > sign_j;
100
101 Ystar_j_ptr.reserve(theta_n(j));
102 alpha_j_ptr.reserve(theta_n(j));
103 theta_j_ptr.reserve(theta_n(j));
104 sign_j.reserve(theta_n(j));
105 /*
106 Allocate proper length for the empty vectors (Ystar_j_ptr, alpha_j_ptr, theta_j_ptr, sign_j)
107 Then we save these pre-allocated vectors inside the outer layer vectors.
108 */
109 theta_Ystar_ptr.push_back(Ystar_j_ptr);
110 theta_alpha_ptr.push_back(alpha_j_ptr);
111 theta_theta_ptr.push_back(theta_j_ptr);
112 theta_sign.push_back(sign_j);
113 }
114
115 for (unsigned int i = 0; i < I; ++i){
116 vector< double* > Ystar_i_ptr;
117 vector< double* > theta1_i_ptr;
118 vector< double* > theta2_i_ptr;
119
120 Ystar_i_ptr.reserve(alpha_n(i));
121 theta1_i_ptr.reserve(alpha_n(i));
122 theta2_i_ptr.reserve(alpha_n(i));
123 /*
124 Allocate proper length for the empty vectors (Ystar_j_ptr, theta1_i_ptr, theta2_i_ptr)
125 Then we save these pre-allocated vectors inside the outer layer vectors.
126 */
127 alpha_Ystar_ptr.push_back(Ystar_i_ptr);
128 alpha_theta1_ptr.push_back(theta1_i_ptr);
129 alpha_theta2_ptr.push_back(theta2_i_ptr);
130 }
131
132
133 for (unsigned int i = 0; i < N; ++i){
134 unsigned int resp = MD(i,0);
135 unsigned int c1 = MD(i,1);
136 unsigned int c2 = MD(i,2);
137 /*
138 Assign the starting values of Ystar, alpha, theta and sign to the above vectors.
139 This is convenient for accessving wanted objects in later updation.
140 */
141 theta_Ystar_ptr[c1].push_back(&Ystar(i,0));
142 theta_alpha_ptr[c1].push_back(&alpha(MD(i,0)));
143 theta_theta_ptr[c1].push_back(&theta(MD(i,2)));
144 theta_sign[c1].push_back(1.0);
145
146 theta_Ystar_ptr[c2].push_back(&Ystar(i,0));
147 theta_alpha_ptr[c2].push_back(&alpha(MD(i,0)));
148 theta_theta_ptr[c2].push_back(&theta(MD(i,1)));
149 theta_sign[c2].push_back(-1.0);
150
151 alpha_Ystar_ptr[resp].push_back(&Ystar(i,0));
152 alpha_theta1_ptr[resp].push_back(&theta(MD(i,1)));
153 alpha_theta2_ptr[resp].push_back(&theta(MD(i,2)));
154 }
155
156
157
158
159
160 // storage matrices (col major order)
161 Matrix<> theta_store;
162 Matrix<> alpha_store;
163 if (storetheta)
164 theta_store = Matrix<>(nsamp, J);
165
166 if (storealpha)
167 alpha_store = Matrix<>(nsamp, I);
168
169
170
171 unsigned int count = 0;
172 // MCMC sampling occurs in this for loop
173 for (unsigned int iter = 0; iter < tot_iter; ++iter){
174
175 // The following three functions are defined in MCMCfcds.h
176 // sample Ystar
177 paircompare_Ystar_update(Ystar, MD, theta, alpha, stream);
178
179 // sample theta
180 paircompare_theta_update(theta, Ystar, MD, alpha, theta_n, theta_eq,
181 theta_ineq,
182 theta_Ystar_ptr,
183 theta_alpha_ptr,
184 theta_theta_ptr,
185 theta_sign,
186 stream);
187
188 if (alphafixed == 0){
189 // sample alpha
190 paircompare_alpha_update(alpha, Ystar, MD, theta, A0, A0a0,
191 alpha_n,
192 alpha_Ystar_ptr,
193 alpha_theta1_ptr,
194 alpha_theta2_ptr,
195 stream);
196 }
197
198 // print results to screen
199 if (verbose > 0 && iter % verbose == 0) {
200 Rprintf("\n\nMCMCpaircompare iteration %i of %i \n", (iter+1), tot_iter);
201 //Rprintf("theta = \n");
202 //for (int j=0; j<J; ++j)
203 //Rprintf("%10.5f\n", theta[j]);
204 }
205
206 // store results
207 if ((iter >= burnin) && ((iter % thin == 0))) {
208
209 // store theta
210 if (storetheta)
211 theta_store(count, _) = theta;
212
213 // store alpha
214 if (storealpha)
215 alpha_store(count, _) = alpha;
216
217 count++;
218 }
219
220 R_CheckUserInterrupt(); // allow user interrupts
221
222 } // end MCMC sampling loop
223
224
225
226
227 // put output back into sampledata
228 Matrix<> output;
229 if(storetheta && ! storealpha) {
230 output = theta_store;
231 } else if (storealpha && ! storetheta){
232 output = alpha_store;
233 } else {
234 output = cbind(theta_store, alpha_store);
235 }
236
237 for (unsigned int i = 0; i < samplesize; ++i)
238 sampledata[i] = output[i];
239
240
241 } // end MCMCpaircompare implementation
242
243
244
245
246
247
248
249
250
251 extern "C" {
252
253 void
cMCMCpaircompare(double * sampledata,const int * samplerow,const int * samplecol,const unsigned int * MDdata,const int * MDrow,const int * MDcol,const int * alphafixed,const int * burnin,const int * mcmc,const int * thin,const int * uselecuyer,const int * seedarray,const int * lecuyerstream,const int * verbose,const double * thetastartdata,const int * thetastartrow,const int * thetastartcol,const double * astartdata,const int * astartrow,const int * astartcol,const double * a0,const double * A0,const double * thetaeqdata,const int * thetaeqrow,const int * thetaeqcol,const double * thetaineqdata,const int * thetaineqrow,const int * thetaineqcol,const int * storealpha,const int * storetheta)254 cMCMCpaircompare(double* sampledata,
255 const int* samplerow,
256 const int* samplecol,
257 const unsigned int* MDdata,
258 const int* MDrow,
259 const int* MDcol,
260 const int* alphafixed,
261 const int* burnin,
262 const int* mcmc,
263 const int* thin,
264 const int *uselecuyer,
265 const int *seedarray,
266 const int *lecuyerstream,
267 const int* verbose,
268 const double* thetastartdata,
269 const int* thetastartrow,
270 const int* thetastartcol,
271 const double* astartdata,
272 const int* astartrow,
273 const int* astartcol,
274 const double* a0,
275 const double* A0,
276 const double* thetaeqdata,
277 const int* thetaeqrow,
278 const int* thetaeqcol,
279 const double* thetaineqdata,
280 const int* thetaineqrow,
281 const int* thetaineqcol,
282 const int* storealpha,
283 const int* storetheta){
284
285
286 // put together matrices
287 const Matrix<unsigned int> MD(*MDrow, *MDcol, MDdata);
288 Matrix<> theta(*thetastartrow, *thetastartcol, thetastartdata);
289 Matrix<> alpha(*astartrow, *astartcol, astartdata);
290 const Matrix<> theta_eq(*thetaeqrow, *thetaeqcol, thetaeqdata);
291 const Matrix<> theta_ineq(*thetaineqrow, *thetaineqcol, thetaineqdata);
292 const int samplesize = (*samplerow) * (*samplecol);
293
294
295 MCMCPACK_PASSRNG2MODEL(MCMCpaircompare_impl, MD, theta, alpha,
296 theta_eq, theta_ineq, *a0, *A0,
297 *alphafixed,
298 *burnin, *mcmc, *thin,
299 *verbose, *storealpha, *storetheta,
300 sampledata, samplesize);
301 }
302 }
303
304
305 #endif
306
307
308
309
310
311
312
313
314
315
316