1 //////////////////////////////////////////////////////////////////////////
2 // cMCMClogit.cc is C++ code to estimate a logistic regression model with
3 // a multivariate normal prior
4 //
5 // Andrew D. Martin
6 // Dept. of Political Science
7 // Washington University in St. Louis
8 // admartin@wustl.edu
9 //
10 // Kevin M. Quinn
11 // Dept. of Government
12 // Harvard University
13 // kevin_quinn@harvard.edu
14 //
15 // This software is distributed under the terms of the GNU GENERAL
16 // PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
17 // file for more information.
18 //
19 // updated to the new version of Scythe 7/25/2004 KQ
20 //
21 // Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
22 // Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
23 // and Jong Hee Park
24 //////////////////////////////////////////////////////////////////////////
25
26
27 #ifndef CMCMCLOGIT_CC
28 #define CMCMCLOGIT_CC
29
30 #include <iostream>
31 #include "MCMCrng.h"
32 #include "MCMCfcds.h"
33 #include "matrix.h"
34 #include "distributions.h"
35 #include "stat.h"
36 #include "la.h"
37 #include "ide.h"
38 #include "smath.h"
39
40 #include <R.h> // needed to use Rprintf()
41 #include <R_ext/Utils.h> // needed to allow user interrupts
42
43 using namespace scythe;
44 using namespace std;
45
46 static double
logit_logpost(const Matrix<> & Y,const Matrix<> & X,const Matrix<> & beta,const Matrix<> & beta_prior_mean,const Matrix<> & beta_prior_prec)47 logit_logpost(const Matrix<>& Y, const Matrix<>& X, const Matrix<>& beta,
48 const Matrix<>& beta_prior_mean,
49 const Matrix<>& beta_prior_prec)
50 {
51 // likelihood
52 const Matrix<> eta = X * beta;
53 const Matrix<> p = 1.0 / (1.0 + exp(-eta));
54 double loglike = 0.0;
55
56 for (unsigned int i = 0; i < Y.rows(); ++ i)
57 loglike += Y(i) * ::log(p(i)) + (1 - Y(i)) * ::log(1 - p(i));
58
59 //prior
60 double logprior = 0.0;
61 if (beta_prior_prec(0) != 0)
62 logprior = lndmvn(beta, beta_prior_mean, invpd(beta_prior_prec));
63
64 return (loglike + logprior);
65 }
66
67 template <typename RNGTYPE>
68 void
MCMClogit_impl(rng<RNGTYPE> & stream,const Matrix<> & Y,const Matrix<> & X,const Matrix<> & tune,Matrix<> & beta,const Matrix<> & b0,const Matrix<> & B0,const Matrix<> & V,unsigned int burnin,unsigned int mcmc,unsigned int thin,unsigned int verbose,Matrix<> & result)69 MCMClogit_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
70 const Matrix<>& X, const Matrix<>& tune, Matrix<>& beta,
71 const Matrix<>& b0, const Matrix<>& B0,
72 const Matrix<>& V, unsigned int burnin, unsigned int mcmc,
73 unsigned int thin, unsigned int verbose,
74 Matrix<>& result)
75 {
76
77 // define constants
78 const unsigned int tot_iter = burnin + mcmc; // total mcmc iterations
79 const unsigned int k = X.cols();
80
81 // proposal parameters
82 const Matrix<> propV = tune * invpd(B0 + invpd(V)) * tune;
83 const Matrix<> propC = cholesky(propV) ;
84
85 double logpost_cur = logit_logpost(Y, X, beta, b0, B0);
86
87 // MCMC loop
88 unsigned int count = 0;
89 unsigned int accepts = 0;
90 for (unsigned int iter = 0; iter < tot_iter; ++iter) {
91
92 // sample beta
93 const Matrix<> beta_can = gaxpy(propC, stream.rnorm(k, 1, 0, 1), beta);
94
95 const double logpost_can = logit_logpost(Y, X, beta_can, b0, B0);
96 const double ratio = ::exp(logpost_can - logpost_cur);
97
98 if (stream.runif() < ratio) {
99 beta = beta_can;
100 logpost_cur = logpost_can;
101 ++accepts;
102 }
103
104 // store values in matrices
105 if (iter >= burnin && ((iter % thin) == 0)) {
106 result(count++, _) = beta;
107 }
108
109 // print output to stdout
110 if(verbose > 0 && iter % verbose == 0){
111 Rprintf("\n\nMCMClogit iteration %i of %i \n", (iter+1), tot_iter);
112 Rprintf("beta = \n");
113 for (unsigned int j = 0; j < k; ++j)
114 Rprintf("%10.5f\n", beta(j));
115 Rprintf("Metropolis acceptance rate for beta = %3.5f\n\n",
116 static_cast<double>(accepts) / static_cast<double>(iter+1));
117 }
118
119 R_CheckUserInterrupt(); // allow user interrupts
120
121 }// end MCMC loop
122 if (verbose > 0){
123 Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
124 Rprintf("The Metropolis acceptance rate for beta was %3.5f",
125 static_cast<double>(accepts) / static_cast<double>(tot_iter));
126 Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
127 }
128 }
129
130 extern "C"{
131
cMCMClogit(double * sampledata,const int * samplerow,const int * samplecol,const double * Ydata,const int * Yrow,const int * Ycol,const double * Xdata,const int * Xrow,const int * Xcol,const int * burnin,const int * mcmc,const int * thin,const double * tunedata,const int * tunerow,const int * tunecol,const int * uselecuyer,const int * seedarray,const int * lecuyerstream,const int * verbose,const double * betastartdata,const int * betastartrow,const int * betastartcol,const double * b0data,const int * b0row,const int * b0col,const double * B0data,const int * B0row,const int * B0col,const double * Vdata,const int * Vrow,const int * Vcol)132 void cMCMClogit(double *sampledata, const int *samplerow,
133 const int *samplecol, const double *Ydata,
134 const int *Yrow, const int *Ycol, const double *Xdata,
135 const int *Xrow, const int *Xcol, const int *burnin,
136 const int *mcmc, const int *thin, const double *tunedata,
137 const int *tunerow, const int *tunecol,
138 const int *uselecuyer, const int *seedarray,
139 const int *lecuyerstream, const int *verbose,
140 const double *betastartdata, const int *betastartrow,
141 const int *betastartcol, const double *b0data,
142 const int *b0row, const int *b0col, const double *B0data,
143 const int *B0row, const int *B0col, const double *Vdata,
144 const int *Vrow, const int *Vcol)
145 {
146
147 // pull together Matrix objects
148 Matrix<> Y(*Yrow, *Ycol, Ydata);
149 Matrix<> X(*Xrow, *Xcol, Xdata);
150 Matrix<> tune(*tunerow, *tunecol, tunedata);
151 Matrix<> beta(*betastartrow, *betastartcol, betastartdata);
152 Matrix<> b0(*b0row, *b0col, b0data);
153 Matrix<> B0(*B0row, *B0col, B0data);
154 Matrix<> V(*Vrow, *Vcol, Vdata);
155
156 Matrix<> result(*samplerow, *samplecol, false);
157 MCMCPACK_PASSRNG2MODEL(MCMClogit_impl, Y, X, tune, beta, b0, B0, V,
158 *burnin, *mcmc, *thin, *verbose, result);
159
160 unsigned int size = *samplecol * *samplerow;
161 for (unsigned int i = 0; i < size; ++i)
162 sampledata[i] = result(i);
163 }
164 }
165
166 #endif
167