1 // cMCMCquantreg.cc is a function that draws from the posterior
2 // distribution of a linear regression model with errors having tauth quantile equal to zero.
3 //
4 // The initial version of this file was generated by the
5 // auto.Scythe.call() function in the MCMCpack R package
6 // written by:
7 //
8 // Andrew D. Martin
9 // Dept. of Political Science
10 // Washington University in St. Louis
11 // admartin@wustl.edu
12 //
13 // Kevin M. Quinn
14 // Dept. of Government
15 // Harvard University
16 // kevin_quinn@harvard.edu
17 //
18 // This software is distributed under the terms of the GNU GENERAL
19 // PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
20 // file for more information.
21 //
22 // Copyright (C) 2009 Andrew D. Martin and Kevin M. Quinn
23 //
24 // This file was initially generated on Wed Apr 1 11:39:12 2009
25 //
26 // The function was rewritten by:
27 //
28 // Craig Reed
29 // Department of Mathematical Sciences
30 // Brunel University
31 // craig.reed@brunel.ac.uk
32 //
33 // CR 9/4/09 [Rewritten function using templates]
34 //
35 // CR 7/12/09 [Placed MCMCmedreg within MCMCquantreg]
36
37 #ifndef CMCMCQUANTREG_CC
38 #define CMCMCQUANTREG_CC
39
40 #include "MCMCrng.h"
41 #include "MCMCfcds.h"
42 #include "matrix.h"
43 #include "distributions.h"
44 #include "stat.h"
45 #include "la.h"
46 #include "ide.h"
47 #include "smath.h"
48 #include "rng.h"
49
50 #include <R.h> // needed to use Rprintf()
51 #include <R_ext/Utils.h> // needed to allow user interrupts
52
53 using namespace std;
54 using namespace scythe;
55
56 /* MCMCquantreg implementation. Takes Matrix<> reference which it
57 * fills with the posterior.
58 */
59 template <typename RNGTYPE>
MCMCquantreg_impl(rng<RNGTYPE> & stream,double tau,Matrix<> & Y,const Matrix<> & X,Matrix<> & beta,const Matrix<> & b0,const Matrix<> & B0,unsigned int burnin,unsigned int mcmc,unsigned int thin,unsigned int verbose,Matrix<> & result)60 void MCMCquantreg_impl (rng<RNGTYPE>& stream, double tau, Matrix<>& Y,
61 const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
62 const Matrix<>& B0,
63 unsigned int burnin, unsigned int mcmc, unsigned int thin,
64 unsigned int verbose,
65 Matrix<>& result)
66 {
67 // define constants
68 const unsigned int tot_iter = burnin + mcmc; //total iterations
69 const unsigned int nstore = mcmc / thin; // number of draws to store
70 const unsigned int k = X.cols ();
71
72 // storage matrices
73 Matrix<> betamatrix (k, nstore);
74
75 // Gibbs sampler
76 unsigned int count = 0;
77 for (unsigned int iter = 0; iter < tot_iter; ++iter) {
78 Matrix<> e = gaxpy(X, (-1*beta), Y);
79 Matrix<> abse = fabs(e);
80 Matrix<> weights = ALaplaceIGaussregress_weights_draw (abse, stream);
81 beta = ALaplaceNormregress_beta_draw (tau, X, Y, weights, b0, B0, stream);
82
83 // store draws in storage matrix
84 if (iter >= burnin && (iter % thin == 0)) {
85 betamatrix(_, count) = beta;
86 ++count;
87 }
88
89 // print output to stdout
90 if(verbose > 0 && iter % verbose == 0) {
91 Rprintf("\n\nMCMCquantreg iteration %i of %i \n",
92 (iter+1), tot_iter);
93 Rprintf("beta = \n");
94 for (unsigned int j=0; j<k; ++j)
95 Rprintf("%10.5f\n", beta(j));
96 }
97
98 R_CheckUserInterrupt(); // allow user interrupts
99
100 } // end MCMC loop
101
102 result = t(betamatrix);
103 } // end MCMCquantreg_impl
104
105 extern "C" {
cMCMCquantreg(double * sampledata,const int * samplerow,const int * samplecol,const double * tau,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 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)106 void cMCMCquantreg(double *sampledata, const int *samplerow,
107 const int *samplecol, const double *tau, const double *Ydata, const int *Yrow,
108 const int *Ycol, const double *Xdata, const int *Xrow,
109 const int *Xcol, const int *burnin, const int *mcmc,
110 const int *thin, const int *uselecuyer, const int *seedarray,
111 const int *lecuyerstream, const int *verbose,
112 const double *betastartdata, const int *betastartrow,
113 const int *betastartcol, const double *b0data,
114 const int *b0row, const int *b0col,
115 const double *B0data, const int *B0row,
116 const int *B0col)
117 {
118 // pull together Matrix objects
119 Matrix<> Y(*Yrow, *Ycol, Ydata);
120 Matrix<> X(*Xrow, *Xcol, Xdata);
121 Matrix<> betastart(*betastartrow, *betastartcol, betastartdata);
122 Matrix<> b0(*b0row, *b0col, b0data);
123 Matrix<> B0(*B0row, *B0col, B0data);
124
125 Matrix<> storagematrix;
126 MCMCPACK_PASSRNG2MODEL(MCMCquantreg_impl, *tau, Y, X, betastart, b0, B0,
127 *burnin, *mcmc, *thin, *verbose,
128 storagematrix);
129
130 const unsigned int size = *samplerow * *samplecol;
131 for (unsigned int i = 0; i < size; ++i)
132 sampledata[i] = storagematrix(i);
133 }
134 }
135
136 #endif
137