mcmcirtkdrob.cc

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

	  z < logfun(R, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
		     Lambda_prior_prec, Lambda_ineq, theta_ineq, 
		     k0, k1, c0, d0, c1, d1, rowindex, colindex))){
    double V = stream.runif();
    if (V < 0.5){
      L = L - (R - L);
    }
    else {
      R = R + (R - L);
    }
    --K;
  }      
  
}


// Radford Neal's (2000) stepping out procedure coded for a logdensity
template<typename RNGTYPE>
static void StepOut(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& m, 
		    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 StepOut().\n");
    exit(1);    
  }
  

  L = x0 - w*U;
  R = L + w;
  const double V = stream.runif();
  int J = static_cast<int>(m*V);
  int K = (m-1) - J; 
  while (J > 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))){
    L = L - w;
    J = J - 1;	   
  }
  while (K > 0 && 
	 (z < logfun(R, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
		     Lambda_prior_prec, Lambda_ineq, theta_ineq, 
		     k0, k1, c0, d0, c1, d1, rowindex, colindex))){
    R = R + w;
    K = K - 1;	   
  }
  
}




// Radford Neal's (2000) Accept procedure coded for a logdensity  
static const bool Accept(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 double& x0, 
			 const double& x1, 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 && x1 >= M) || (x0 >= M && x1 < M)){
      D = true;
    }
    if (x1 < M){
      Rhat = M;
    }
    else {
      Lhat = M;
    }
    
    if (D && z >= logfun(Lhat, X, Lambda, theta, delta0, delta1, 
			 Lambda_prior_mean, Lambda_prior_prec, 
			 Lambda_ineq, theta_ineq, k0, k1, c0, d0, 
			 c1, d1, rowindex, colindex) && 
	z >=  logfun(Rhat, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
		     Lambda_prior_prec, Lambda_ineq, theta_ineq, 
		     k0, k1, c0, d0, c1, d1, rowindex, colindex)){
      return(false);
    }    
  }
  return(true);
}


// Radford Neal's (2000) shrinkage procedure coded for a log density
// assumes the doubling procedure has been used to find L and R
template <typename RNGTYPE>
static double shrinkageDoubling(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, 
				rng<RNGTYPE>& stream, const double& L, 
				const double& R,
				const int& param){
  
  double Lbar = L;
  double Rbar = R;
  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 shrinkageDoubling().\n");
    exit(1);    
  }

  for (;;){
    const double U = stream.runif();
    const double x1 = Lbar + U*(Rbar - Lbar);
    if (z <= logfun(x1, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
		    Lambda_prior_prec, Lambda_ineq, theta_ineq, 
		    k0, k1, c0, d0, c1, d1, rowindex, colindex) &&
	Accept(logfun, X, Lambda, theta, delta0, delta1, Lambda_prior_mean, 
	       Lambda_prior_prec, Lambda_ineq, theta_ineq, 
	       k0, k1, c0, d0, c1, d1, rowindex, colindex, z, w, 
	       x0, x1, L, R)){
      return(x1);
    }
    if (x1 < x0){
      Lbar = x1;
    }
    else {
      Rbar = x1;
    }
  } // end infinite loop
}



// Radford Neal's (2000) shrinkage procedure coded for a log density
// assumes the stepping out procedure has been used to find L and R
template <typename RNGTYPE>
static double shrinkageStep(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, 
			    rng<RNGTYPE>& stream, const double& L, 
			    const double& R,
			    const int& param){
  
  double Lbar = L;
  double Rbar = R;
  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 shrinkageDoubling().\n");
    exit(1);    
  }

  for (;;){
    const double U = stream.runif();
    const double x1 = Lbar + U*(Rbar - Lbar);
    if (z <= logfun(x1, X, Lambda, theta, delta0, delta1, Lambda_prior_mean,
		    Lambda_prior_prec, Lambda_ineq, theta_ineq, 
		    k0, k1, c0, d0, c1, d1, rowindex, colindex) ){
      return(x1);
    }
    if (x1 < x0){
      Lbar = x1;
    }
    else {
      Rbar = x1;
    }
  } // end infinite loop
}













template <typename RNGTYPE>
void MCMCirtKdRob_impl(rng<RNGTYPE>& stream,
		       const Matrix<int>& X, 
		       Matrix<>& Lambda,
		       Matrix<>& theta,
		       const Matrix<>& Lambda_eq, const Matrix<>& Lambda_ineq,
		       const Matrix<>& theta_eq, const Matrix<>& theta_ineq,
		       const Matrix<>& Lambda_prior_mean,
		       const Matrix<>& Lambda_prior_prec,
		       const int* burnin, const int* mcmc,  const int* thin,
		       const int* verbose, 
		       const int* method_step,
		       const double* theta_w, const int* theta_p,
		       const double* lambda_w, const int* lambda_p,
		       const double* delta0_w, const int* delta0_p,
		       const double* delta1_w, const int* delta1_p, 
		       const double * delta0start, const double* delta1start,
		       const double* k0, const double* k1,
		       const double* c0, const double* c1,
		       const double* d0, const double* d1,
		       const int* storeitem,
		       const int* storeability,
		       double* sampledata, const int* samplerow, 
		       const int* samplecol
		       ){



    // constants 
    const int K = X.cols();  // number of items
    const int N = X.rows();  // number of subjects
    const int D = Lambda.cols();  // number of dimensions + 1
    const int tot_iter = *burnin + *mcmc;  
    const int nsamp = *mcmc / *thin;
    // const Matrix<double> Lambda_free_indic = Matrix<double>(K, D);
    //for (int i=0; i<(K*D); ++i){
    //  if (Lambda_eq[i] == -999) Lambda_free_indic[i] = 1.0;
    //}

    // starting values
    //  Matrix<double> theta = Matrix<double>(N, D-1);
    Matrix<> onesvec = ones<double>(N, 1);
    onesvec = onesvec * -1.0;