mcmcirtkdrob.cc

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

// MCMCirtKdRob.cc is C++ code to estimate a robust K-dimensional
// item response model 
//
// 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
// 
// 2/18/2005 KQ
// 8/1/2007  ported to Scythe 1.0.2 KQ
//


#ifndef MCMCIRTKDROB_CC
#define MCMCIRTKDROB_CC

#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 "MCMCmnl.h"

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


typedef Matrix<double,Row,View> rmview;

using namespace std;
using namespace scythe;


/* Equal probability sampling; without-replacement case */
// pulled from R-2.0.1/src/main/random.c lines 352-364
// slightly modified by KQ 2/21/2005
// x: n array of original indices running from 0 to (n-1)
// y: k array of samled indices
// k: length of y (must be <= n)
// n: length of x (must be >= k)
template <typename RNGTYPE>
static void SampleNoReplace(const int k, int n, int *y, int *x, 
			    rng<RNGTYPE>& stream){
  
  for (int i = 0; i < n; i++)
    x[i] = i;
  for (int i = 0; i < k; i++) {
    int j = static_cast<int>(n * stream.runif());
    y[i] = x[j];
    x[j] = x[--n];
  }
}



// full conditional distribution for a single element of Lambda
// this single element is Lambda(rowindex, colindex)
static double Lambda_logfcd(const double& lam_ij,
			    const Matrix<int>& X,
			    const Matrix<>& Lambda,
			    const Matrix<>& theta, 
			    const double& delta0,
			    const double& delta1, 
			    const Matrix<>& Lambda_prior_mean, 
			    const Matrix<>& Lambda_prior_prec,
			    const Matrix<>& Lambda_ineq,
			    const Matrix<>& theta_ineq,
			    const double& k0,
			    const double& k1,
			    const double& c0,
			    const double& d0,
			    const double& c1,
			    const double& d1,
			    const int& rowindex,
			    const int& colindex){


  const int D = Lambda.cols();
  
  // check to see if inequality constraint is satisfied and 
  // evaluate prior
    
  double logprior = 0.0;  
  if (Lambda_ineq(rowindex,colindex) * lam_ij < 0){
    return log(0.0);
  } 
  if (Lambda_prior_prec(rowindex,colindex) != 0){
    logprior += lndnorm(lam_ij, 
			Lambda_prior_mean(rowindex,colindex), 
			sqrt(1.0 / Lambda_prior_prec(rowindex,colindex)));
  } 
  
  // prior is uniform on hypersphere with radius 10
  /*
    if (Lambda_ineq(rowindex,colindex) * lam_ij < 0){
    return log(0.0);
    } 
    double absqdist = 0.0;
    for (int i=0; i<D; ++i){
    if (i==colindex){
    absqdist += ::pow(lam_ij, 2);
    }
    else{
    absqdist += ::pow(Lambda(rowindex,i), 2);
    }
    }
    if (absqdist > 100.0){
    return log(0.0);
    }
    const double logprior = 0.0;
  */

  // likelihood
  double loglike = 0.0;
  const int N = X.rows();
  for (int i=0; i<N; ++i){
    if (X(i,rowindex) != -999){
      double eta = 0.0;
      for (int j=0; j<D; ++j){
	if (j==colindex){
	  eta += theta(i,j) * lam_ij;
	}
	else{
	  eta += theta(i,j) * Lambda(rowindex,j);
	}
      } 
      const double p = delta0 + (1 - delta0 - delta1) * 
	(1.0 / (1.0 + exp(-1*eta))); 
      loglike += X(i,rowindex) * log(p) + (1.0 - X(i,rowindex)) * log(1.0 - p);
    } 
  }
  
  return (loglike + logprior);
}



// full conditional for a single element of theta
// this single element is theta(rowindex, colindex)
static double theta_logfcd(const double& t_ij,
			   const Matrix<int>& X, 
			   const Matrix<>& Lambda,
			   const Matrix<>& theta, 
			   const double& delta0,
			   const double& delta1, 
			   const Matrix<>& Lambda_prior_mean, 
			   const Matrix<>& Lambda_prior_prec,
			   const Matrix<>& Lambda_ineq,
			   const Matrix<>& theta_ineq,
			   const double& k0,
			   const double& k1,
			   const double& c0,
			   const double& d0,
			   const double& c1,
			   const double& d1,
			   const int& rowindex,
			   const int& colindex){
 
  const int D = Lambda.cols();     
  // evaluate prior  
  if (theta_ineq(rowindex,colindex-1) * t_ij < 0){
    return log(0.0);
  } 
  const double logprior = lndnorm(t_ij, 0.0, 1.0); 
  
  // prior is uniform on unit circle
  /*
    if (theta_ineq(rowindex,colindex-1) * t_ij < 0){
    return log(0.0);
    } 
    double tsqdist = 0.0;
    for (int i=0; i<(D-1); ++i){
    if (i==(colindex-1)){
    tsqdist += ::pow(t_ij, 2);
    }
    else{
    tsqdist += ::pow(theta(rowindex,(i-1)), 2);
    }
    }
    if (tsqdist > 1.0){
    return log(0.0);
    }
    const double logprior = 1.0;
  */

  // likelihood
  double loglike = 0.0;
  const int K = X.cols();
  for (int i=0; i<K; ++i){
    if (X(rowindex,i) != -999){
      double eta = 0.0;
      for (int j=0; j<D; ++j){
	if (j==colindex){
	  eta += t_ij * Lambda(i,j);
	}
	else{
	  eta += theta(rowindex,j) * Lambda(i,j);
	}
      }       
      const double p = delta0 + (1 - delta0 - delta1) * 
	(1.0 / (1.0 + exp(-1*eta)));       
      loglike += X(rowindex,i) * log(p) + (1.0 - X(rowindex,i)) * log(1.0 - p);
    } 
  }
  
  return (loglike + logprior);
}



// full conditional for delta0
static double delta0_logfcd(const double& delta0,
			    const Matrix<int>& X, 
			    const Matrix<>& Lambda,
			    const Matrix<>& theta, 
			    const double& junk,
			    const double& delta1, 
			    const Matrix<>& Lambda_prior_mean, 
			    const Matrix<>& Lambda_prior_prec,
			    const Matrix<>& Lambda_ineq,
			    const Matrix<>& theta_ineq,
			    const double& k0,
			    const double& k1,
			    const double& c0,
			    const double& d0,
			    const double& c1,
			    const double& d1,
			    const int& rowindex,
			    const int& colindex){
  
  // evaluate prior
  if (delta0 >=k0 || delta0 <=0){
    return log(0.0);
  }
  const double logprior = lndbeta1(delta0 * (1.0/k0), c0, d0); 
  
  // likelihood
  double loglike = 0.0;
  const int D = Lambda.cols();
  const int N = X.rows();
  const int K = X.cols();
  for (int i=0; i<N; ++i){
    for (int k=0; k<K; ++k){
      if (X(i,k) != -999){
	double eta = 0.0;
	for (int j=0; j<D; ++j){
	  eta += theta(i,j) * Lambda(k,j);
	} 
	const double p = delta0 + (1 - delta0 - delta1) * 
	  (1.0 / (1.0 + exp(-1*eta))); 
	loglike += X(i,k) * log(p) + (1.0 - X(i,k)) * log(1.0 - p);
      }
    } 
  }
  
  return (loglike + logprior);
}



// full conditional for delta1
static double delta1_logfcd(const double& delta1, 
			    const Matrix<int>& X, 
			    const Matrix<>& Lambda,
			    const Matrix<>& theta, 
			    const double& delta0,
			    const double& junk, 
			    const Matrix<>& Lambda_prior_mean, 
			    const Matrix<>& Lambda_prior_prec,
			    const Matrix<>& Lambda_ineq,
			    const Matrix<>& theta_ineq,
			    const double& k0,
			    const double& k1,
			    const double& c0,
			    const double& d0,
			    const double& c1,
			    const double& d1,
			    const int& rowindex,
			    const int& colindex){
  
  // evaluate prior
  if (delta1 >=k1 || delta1 <=0){
    return log(0.0);
  }
  const double logprior = lndbeta1(delta1 * (1.0/k1), c1, d1); 
  
  // likelihood
  double loglike = 0.0;
  const int D = Lambda.cols();
  const int N = X.rows();
  const int K = X.cols();
  for (int i=0; i<N; ++i){
    for (int k=0; k<K; ++k){
      if (X(i,k) != -999){
	double eta = 0.0;
	for (int j=0; j<D; ++j){
	  eta += theta(i,j) * Lambda(k,j);
	} 
	const double p = delta0 + (1 - delta0 - delta1) * 
	  (1.0 / (1.0 + exp(-1*eta))); 
	loglike += X(i,k) * log(p) + (1.0 - X(i,k)) * log(1.0 - p);
      }
    } 
  }
  
  return (loglike + logprior);
}





// Radford Neal's (2000) doubling procedure coded for a logdensity
template<typename RNGTYPE>
static void doubling(double (*logfun)(const double&,
				      const Matrix<int>&,
				      const Matrix<>&,
				      const Matrix<>&,
				      const double&, 
				      const double&,
				      const Matrix<>&,
				      const Matrix<>&,
				      const Matrix<>&,
				      const Matrix<>&,
				      const double&, 
				      const double&,
				      const double&, 
				      const double&,
				      const double&, 
				      const double&,
				      const int&,
				      const int&), 
		     const Matrix<int>& X,
		     const Matrix<>& Lambda,
		     const Matrix<>& theta,
		     const double& delta0,
		     const double& delta1, 
		     const Matrix<>& Lambda_prior_mean, 
		     const Matrix<>& Lambda_prior_prec,
		     const Matrix<>& Lambda_ineq,
		     const Matrix<>& theta_ineq,
		     const double& k0,
		     const double& k1,
		     const double& c0,
		     const double& d0,
		     const double& c1, 
		     const double& d1,
		     const int& rowindex,
		     const int& colindex,
		     const double& z, 
		     const double& w, const int& p, 
		     rng<RNGTYPE>& stream, double& L, double& R, 
		     const int& param){
  
  const double U = stream.runif();
  double x0;
  if (param==0){ // Lambda
    x0 = Lambda(rowindex, colindex);
  }
  else if (param==1){ // theta
    x0 = theta(rowindex, colindex);
  }
  else if (param==2){ // delta0
    x0 = delta0;
  }
  else if (param==3){ // delta1
    x0 = delta1;
  }
  else {
    Rprintf("\nERROR: param not in {0,1,2,3} in doubling().\n");
    exit(1);    
  }

  L = x0 - w*U;
  R = L + w;
  int K = p;
  while (K > 0 && 
	 (z < logfun(L, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
		     Lambda_prior_prec, Lambda_ineq, theta_ineq, 
		     k0, k1, c0, d0, c1, d1, rowindex, colindex) |