Logo Search packages:      
Sourcecode: mcmcpack version File versions

MCMCregress.cc

// MCMCregress.cc is a program that simualates draws from the posterior
// density of a linear regression model with Gaussian errors.
//
// 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
// admartin@wustl.edu
//
// Kevin M. Quinn
// Dept. of Government
// Harvard University
// kevin_quinn@harvard.edu
// 
// 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 Fri Jul 23 15:07:21 2004
//
// ADM and KQ 10/10/2002 [ported to Scythe0.3]
// ADM 6/2/04 [re-written using template]
// KQ 6/18/04 [modified to meet new developer specification]
// ADM 7/22/04 [modified to work with new Scythe and rngs]

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

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

using namespace SCYTHE;
using namespace std;

extern "C" {

   // simulate from posterior density and return an mcmc by parameters
   // matrix of the posterior density sample
   void MCMCregress(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 int *lecuyer, 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
     Matrix <double> Y = r2scythe(*Yrow, *Ycol, Ydata);
     Matrix <double> X = r2scythe(*Xrow, *Xcol, Xdata);
     Matrix <double> betastart = r2scythe(*betastartrow,
       *betastartcol, betastartdata);
     Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
     Matrix <double> B0 = r2scythe(*B0row, *B0col, B0data);

     // define constants and form cross-product matrices
     const int tot_iter = *burnin + *mcmc; // total number of mcmc iterations
     const int nstore = *mcmc / *thin; // number of draws to store
     const int k = X.cols ();
     const Matrix <double> XpX = crossprod(X);
     const Matrix <double> XpY = t(X) * Y;

     // storage matrices
     Matrix <double> betamatrix (k, nstore);
     Matrix <double> sigmamatrix (1, nstore);

     // initialize rng stream
     rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);

     // set starting values
     Matrix <double> beta = betastart;

     // Gibbs sampler
     int count = 0;
     for (int iter = 0; iter < tot_iter; ++iter) {
       double sigma2 = NormIGregress_sigma2_draw (X, Y, beta, *c0,
         *d0, stream);
       beta = NormNormregress_beta_draw (XpX, XpY, b0, B0, sigma2,
         stream);  
         
       // store draws in storage matrix (or matrices)
       if (iter >= *burnin && (iter % *thin == 0)) {
         sigmamatrix (0, count) = sigma2;
         for (int j = 0; j < k; j++)
           betamatrix (j, count) = beta[j];
         ++count;
       }
       
       // print output to stdout
       if(*verbose > 0 && iter % *verbose == 0) {
         Rprintf("\n\nMCMCregress iteration %i of %i \n",
           (iter+1), tot_iter);
         Rprintf("beta = \n");
         for (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

     delete stream; // clean up random number stream

     // load draws into sample array
     Matrix <double> storeagematrix = cbind (t (betamatrix), t (sigmamatrix));     
     const int size = *samplerow * *samplecol;
     for(int i = 0; i < size; ++i)
       sampledata[i] = storeagematrix[i];

   } // end MCMCregress 
} // end extern "C"

Generated by  Doxygen 1.6.0   Back to index