Logo Search packages:      
Sourcecode: mcmcpack version File versions


// 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
// 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 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]

#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" {

   // 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 *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, 
     Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
     Matrix <double> B0 = r2scythe(*B0row, *B0col, B0data);
     Matrix <double> t_X = t(X);
     // define constants
     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 int N = X.rows ();
     const Matrix <double> XpX = crossprod(X);
     // storage matrix or 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;
     int count = 0;
     Matrix <double> Z = Y;
     for(int iter = 0; iter < tot_iter; ++iter){
       double sigma2 = NormIGregress_sigma2_draw (X, Z, beta, *c0, *d0, 
       Matrix <double> Z_mean = X * beta;
       for (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 <double> 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 % 1 == 0)) {

         sigmamatrix (0, count) = sigma2;
         for (int j = 0; j < k; j++)
           betamatrix (j, count) = beta[j];
       // 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 (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 MCMCtobit 
} // end extern "C"

Generated by  Doxygen 1.6.0   Back to index