mcmctobit.cc

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

// MCMCtobit.cc is a program that simualates draws from the posterior
// density of a linear regression model with Gaussian errors when the 
// dependent variable is censored from below and/or above.
//
// The initial version of this file was generated by the
// auto.Scythe.call() function in the MCMCpack R package
// written by:
//
// Andrew D. Martin
// Dept. of Political Science
// Washington University in St. Louis
// [e-mail]
//
// Kevin M. Quinn
// Dept. of Government
// Harvard University
// [e-mail]
// 
// This software is distributed under the terms of the GNU GENERAL
// PUBLIC LICENSE Version 2, June 1991.  See the package LICENSE
// file for more information.
//
// Copyright (C) 2004 Andrew D. Martin and Kevin M. Quinn
// 
// This file was initially generated on Tue Sep 14 00:50:08 2004
// ADM and KQ 10/10/2002 [ported to Scythe0.3]
// BG 09/18/2004 [updated to new specification, added above censoring]
// ADM 7/7/2007 [updated to Scythe 1.0.X]

#ifndef MCMCTOBIT_CC
#define MCMCTOBIT_CC

#include "MCMCrng.h"
#include "MCMCfcds.h"
#include "matrix.h"
#include "distributions.h"
#include "stat.h"
#include "la.h"
#include "ide.h"
#include "smath.h"
#include "rng.h"
#include "mersenne.h"
#include "lecuyer.h"

#include <R.h>           // needed to use Rprintf()
#include <R_ext/Utils.h> // needed to allow user interrupts

using namespace std;
using namespace scythe;

/* MCMCtobit implemenation.  Takes Matrix<> reference which it
 * fills with the posterior. 
 */
template <typename RNGTYPE>
void MCMCtobit_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
    const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
    const Matrix<>& B0, double c0, double d0, double below, double above,
    unsigned int burnin, unsigned int mcmc, unsigned int thin, 
    unsigned int verbose, Matrix<>& result) {
     
     // define constants
     const unsigned int tot_iter = burnin + mcmc;  // total number of mcmc iterations
     const unsigned int nstore = mcmc / thin;      // number of draws to store
     const unsigned int k = X.cols();
     const unsigned int N = X.rows();
     const Matrix <> XpX = crossprod(X);
     
     // storage matrix or matrices
     Matrix <> betamatrix (k, nstore);
     Matrix <> sigmamatrix (1, nstore);
     
     ///// MCMC SAMPLING OCCURS IN THIS FOR LOOP
     int count = 0;
     Matrix <> Z = Y;
     for(unsigned int iter = 0; iter < tot_iter; ++iter){
        double sigma2 = NormIGregress_sigma2_draw (X, Z, beta, c0, d0, 
						  stream);
        Matrix <> Z_mean = X * beta;
       
        for (unsigned int i=0; i<N; ++i) {
           if (Y[i] <= below)
	          Z[i] = stream.rtanorm_combo(Z_mean[i], sigma2, below);
	       if (Y[i] >= above)
	          Z[i] = stream.rtbnorm_combo(Z_mean[i], sigma2, above);
        }
        Matrix <> XpZ = t(X) * Z;
        beta = NormNormregress_beta_draw (XpX, XpZ, b0, B0, sigma2, stream);  
       
        // store draws in storage matrix (or matrices)
        if (iter >= burnin && (iter % thin == 0)) {
          sigmamatrix (0, count) = sigma2;
          betamatrix(_, count) = beta;
          ++count;
        }
       
        // print output to stdout
        if(verbose > 0 && iter % verbose == 0) {
           Rprintf("\n\nMCMCtobit iteration %i of %i \n",
		      (iter+1), tot_iter);
           Rprintf("beta = \n");
           for (unsigned int j=0; j<k; ++j)
              Rprintf("%10.5f\n", beta[j]);
           Rprintf("sigma2 = %10.5f\n", sigma2);
        }

        R_CheckUserInterrupt(); // allow user interrupts  
     } // end MCMC loop
     result = cbind (t(betamatrix), t(sigmamatrix));
}

extern "C" {

   // MCMCtobit is a linear regression model with a censored dependent variable
   void MCMCtobit(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 double *below, const double *above, 
		  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, const double *c0, 
		  const double *d0) {
   
     // pull together Matrix objects
     const Matrix <> Y(*Yrow, *Ycol, Ydata);
     const Matrix <> X(*Xrow, *Xcol, Xdata);
     Matrix <double> betastart(*betastartrow, *betastartcol, 
					  betastartdata);
     const Matrix <> b0(*b0row, *b0col, b0data);
     const Matrix <> B0(*B0row, *B0col, B0data);
         
     Matrix<> storagematrix;
     MCMCPACK_PASSRNG2MODEL(MCMCtobit_impl, Y, X, betastart, b0, B0,
       *c0, *d0, *below, *above, *burnin, *mcmc, *thin, *verbose,
       storagematrix);

    const unsigned int size = *samplerow * *samplecol;
    for (unsigned int i=0; i<size; ++i)
      sampledata[i] = storagematrix(i);
    }
}

#endif