Logo Search packages:      
Sourcecode: mcmcpack version File versions  Download package

MCMCmnlslice.cc

// MCMCmnlslice.cc DESCRIPTION HERE
//
// 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 Wed Dec 29 15:27:40 2004
// REVISION HISTORY

#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;

static double mnl_logpost(const Matrix<double>& Y, const Matrix<double>& X, 
                    const Matrix<double>& beta,
                    const Matrix<double>& beta_prior_mean, 
                    const Matrix<double>& beta_prior_prec){

  //  likelihood
  double loglike = 0.0;
  Matrix<double> numer = exp(X * beta);
  numer = reshape(numer, Y.rows(), Y.cols());
  double *denom = new double[Y.rows()];
  for (int i=0; i<Y.rows(); ++i){
    denom[i] = 0.0;
    for (int j=0; j<Y.cols(); ++j){
      if (Y(i,j) != -999){
      denom[i] += numer(i,j);
      }
    }
    for (int j=0; j<Y.cols(); ++j){
      if (Y(i,j) == 1.0){
      loglike += std::log(numer(i,j) / denom[i]);
      }
    }
  }
  
  delete [] denom;

  // prior
  double logprior = 0.0;
  if (beta_prior_prec(0,0) != 0){
    logprior = lndmvn(beta, beta_prior_mean, invpd(beta_prior_prec));
  }

  return (loglike + logprior);
}




// eventually all of the slice sampling functions should be made more 
// general and put in MCMCfcds.{h cc}
//
// Radford Neal's (2000) doubling procedure coded for a logdensity
  static void doubling(double (*logfun)(const Matrix<double>&,
                              const Matrix<double>&,
                              const Matrix<double>&,
                              const Matrix<double>&,
                              const Matrix<double>&), 
                   const Matrix<double>& beta, const int& index, 
                   const double& z, 
                   const double& w, const int& p, 
                   const Matrix<double>& Y,
                   const Matrix<double>& X,
                   const Matrix<double>& beta_prior_mean,
                   const Matrix<double>& beta_prior_prec,
                   rng* stream, double& L, double& R){
    
  const double U = stream->runif();
  const double x0 = beta[index];
  Matrix<double> beta_L = beta;
  Matrix<double> beta_R = beta;
  L = x0 - w*U;
  beta_L[index] = L;
  R = L + w;
  beta_R[index] = R;
  int K = p;
  while (K > 0 && 
       (z < logfun(Y, X, beta_L, beta_prior_mean, beta_prior_prec) ||
        z < logfun(Y, X, beta_R, beta_prior_mean, beta_prior_prec))){
    double V = stream->runif();
    if (V < 0.5){
      L = L - (R - L);
      beta_L[index] = L;
    }
    else {
      R = R + (R - L);
      beta_R[index] = R;
    }
    --K;
  }  
}

// Radford Neal's (2000) Accept procedure coded for a logdensity
static const bool Accept(double (*logfun)(const Matrix<double>&,
                                const Matrix<double>&,
                                const Matrix<double>&,
                                const Matrix<double>&,
                                const Matrix<double>&), 
                   const Matrix<double>& beta, const int& index, 
                   const double x0, 
                   const double& z, const double& w, 
                   const Matrix<double>& Y,
                   const Matrix<double>& X,
                   const Matrix<double>& beta_prior_mean,
                   const Matrix<double>& beta_prior_prec,
                   const double& L, const double& R){
  
  double Lhat = L;
  double Rhat = R;
  bool D = false;
  while ((Rhat - Lhat ) > 1.1 * w){
    double M = (Lhat + Rhat) / 2.0;
    if ( (x0 < M && beta[index] >= M) || (x0 >= M && beta[index] < M)){
      D = true;
    }
    if (beta[index] < M){
      Rhat = M;
    }
    else {
      Lhat = M;
    }

    Matrix<double> beta_L = beta;
    Matrix<double> beta_R = beta;
    beta_L[index] = Lhat;
    beta_R[index] = Rhat;

    if (D && z >= logfun(Y, X, beta_L, beta_prior_mean, beta_prior_prec) && 
      z >=  logfun(Y, X, beta_R, beta_prior_mean, beta_prior_prec)){
      return(false);
    }    
  }
  return(true);
}


// Radford Neal's (2000) shrinkage procedure coded for a log density
static double shrinkage(double (*logfun)(const Matrix<double>&,
                               const Matrix<double>&,
                               const Matrix<double>&,
                               const Matrix<double>&,
                               const Matrix<double>&),
                  const Matrix<double>& beta, const int& index, 
                  const double& z, const double& w, 
                  const Matrix<double>& Y,
                  const Matrix<double>& X,
                  const Matrix<double>& beta_prior_mean,
                  const Matrix<double>& beta_prior_prec,
                  rng* stream, const double& L, const double& R){
  
  double Lbar = L;
  double Rbar = R;
  Matrix<double> beta_x1 = beta;
  const double x0 = beta[index]; 
  for (;;){
    const double U = stream->runif();
    const double x1 = Lbar + U*(Rbar - Lbar);
    beta_x1[index] = x1;
    if (z < logfun(Y, X, beta_x1, beta_prior_mean, beta_prior_prec) &&
      Accept(logfun, beta_x1, index, x0, z, w, 
             Y, X, beta_prior_mean, beta_prior_prec, L, R)){
      return(x1);
    }
    if (x1 < x0){
      Lbar = x1;
    }
    else {
      Rbar = x1;
    }
  } // end infinite loop
}






extern "C" {

   // MCMC sampling for MNL model via slice sampling
   void MCMCmnlslice(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 *Vdata, const int *Vrow, 
                 const int *Vcol) {
     
     // pull together Matrix objects
     // REMEMBER TO ACCESS PASSED ints AND doubles PROPERLY
     const Matrix <double> Y = r2scythe(*Yrow, *Ycol, Ydata);
     const Matrix <double> X = r2scythe(*Xrow, *Xcol, Xdata);
     Matrix <double> beta = r2scythe(*betastartrow, *betastartcol, 
                             betastartdata);     
     const Matrix <double> b0 = r2scythe(*b0row, *b0col, b0data);
     const Matrix <double> B0 = r2scythe(*B0row, *B0col, B0data);
     const Matrix <double> V = r2scythe(*Vrow, *Vcol, Vdata);
     
     // 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();

     // storage matrix or matrices
     Matrix<double> storemat(nstore, k);

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

     // proposal parameters
     const Matrix<double> propV = invpd(B0 + invpd(V));
     const Matrix<double> w_init = ones<double>(k, 1);
     for (int i=0; i<k; ++i)
       w_init[i] = sqrt(propV(i,i))*0.05;

     // starting values
     double L = -1.0;
     double R = 1.0;


     const int warmup_iter = 100;
     const int warmup_burnin = 10;
     const int p_init = 15;
     const Matrix<double> widthmat(warmup_iter - warmup_burnin, k);
     // warm up sampling to choose the slice sampling parameters
     for (int iter=0; iter<warmup_iter; ++iter){
       for (int index=0; index<k; ++index){
      
       double funval = mnl_logpost(Y, X, beta, b0, 
                             B0);
       double z = funval - stream->rexp(1.0);
       doubling(&mnl_logpost, beta, index, z, w_init[index], p_init, Y, X, 
              b0, B0, stream, L, R);
       
       beta[index] = shrinkage(&mnl_logpost, beta, index, z, 
                         w_init[index], Y, X, b0,
                         B0, stream, L, R);
       if (iter >= warmup_burnin){
         widthmat(iter-warmup_burnin, index) =  R - L;
       }
       }
     }
     const Matrix<double> w = meanc(widthmat);
     Matrix<int> p = ones<int>(k,1);
     for (int index=0; index<k; ++index){
       int p_temp = 2;
       while ((w[index] * pow(2.0, p_temp) ) < max(widthmat(_,index))){
       ++p_temp;
       } 
       p[index] = p_temp + 1;       
     }
     


     int count = 0;
     ///// REAL MCMC SAMPLING OCCURS IN THIS FOR LOOP
     for(int iter = 0; iter < tot_iter; ++iter){
       for (int index=0; index<k; ++index){
       
       double funval = mnl_logpost(Y, X, beta, b0, 
                             B0);
       double z = funval - stream->rexp(1.0);
       doubling(&mnl_logpost, beta, index, z, w[index], p[index], 
              Y, X, b0, B0, stream, L, R);  
       beta[index] = shrinkage(&mnl_logpost, beta, index, z, 
                         w[index], Y, X, b0,
                         B0, stream, L, R);
       }


       // store draws in storage matrix (or matrices)
       if(iter >= *burnin && (iter % *thin == 0)){
       for (int j = 0; j < k; j++)
         storemat(count, j) = beta[j];
       ++count;
       }

       // print output to stdout
       if(*verbose > 0 && iter % *verbose == 0){
         Rprintf("\n\nMCMCmnl slice 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]);
       }

       R_CheckUserInterrupt(); // allow user interrupts
       
     } // end MCMC loop

     delete stream; // clean up random number stream

     // load draws into sample array
     const int size = *samplerow * *samplecol;
     for(int i = 0; i < size; ++i)
       sampledata[i] = storemat[i];

   } // end MCMCmnlslice 
} // end extern "C"

Generated by  Doxygen 1.6.0   Back to index