mcmcirtkdrob.cc

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

    theta = cbind(onesvec, theta);
    double delta0 = *delta0start;
    double delta1 = *delta1start;

    // index arrays (used for the random order of the sampling)
    // OLD
    //int K_array[K];
    //int N_array[N];
    //int D_array[D];
    //int Dm1_array[D-1];
    //int K_inds_perm[K];
    //int N_inds_perm[N];
    //int D_inds_perm[D];
    //int Dm1_inds_perm[D-1];
    // NEW
    int* K_array = new int[K];
    int* N_array = new int[N];
    int* D_array = new int[D];
    int* Dm1_array = new int[D-1];
    int* K_inds_perm = new int[K];
    int* N_inds_perm = new int[N];
    int* D_inds_perm = new int[D];
    int* Dm1_inds_perm = new int[D-1];


    // storage matrices (row major order)
    Matrix<> Lambda_store;
    if (storeitem[0]==1){
      Lambda_store = Matrix<double>(nsamp, K*D);
    }
    Matrix<> theta_store;
    if (*storeability==1){
      theta_store = Matrix<double>(nsamp, N*D);
    }
    Matrix<> delta0_store(nsamp, 1);
    Matrix<> delta1_store(nsamp, 1);

    ///////////////////
    // Slice Sampler //
    ///////////////////
    int count = 0;  
    for (int iter=0; iter < tot_iter; ++iter){

      double L, R, w, funval, z;
      int p;

      // sample theta
      int param = 1;    
      SampleNoReplace(N, N, N_inds_perm, N_array, stream);   
      for (int ii=0; ii<N; ++ii){
	int i = N_inds_perm[ii];
	SampleNoReplace(D-1, D-1, Dm1_inds_perm, Dm1_array, stream);    
	for (int jj=0; jj<(D-1); ++jj){
	  int j = Dm1_inds_perm[jj]+1;
	  if (theta_eq(i,j-1) == -999){
	    w = *theta_w;
	    p = *theta_p;
	    L = -1.0;
	    R = 1.0;
	    funval = theta_logfcd(theta(i,j), X, Lambda, 
				  theta, delta0,
				  delta1, Lambda_prior_mean, 
				  Lambda_prior_prec,
				  Lambda_ineq, theta_ineq, 
				  *k0, *k1, *c0, *d0,
				  *c1, *d1, i, j);
	    z = funval - stream.rexp(1.0);

	    if (*method_step == 1){
	      StepOut(&theta_logfcd, X, Lambda, theta, delta0, delta1, 
		      Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
		      theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, i, j, z, 
		      w, p, stream, L, R, param);
	    
	      theta(i,j) = shrinkageStep(&theta_logfcd, X, Lambda, theta,
					 delta0, delta1, Lambda_prior_mean, 
					 Lambda_prior_prec, Lambda_ineq,
					 theta_ineq, *k0, *k1, 
					 *c0, *d0, *c1, 
					 *d1, i, j, z, w, stream, 
					 L, R, param);
	    }
	    else{
	      doubling(&theta_logfcd, X, Lambda, theta, delta0, delta1, 
		       Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
		       theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, i, j, z, 
		       w, p, stream, L, R, param);
	      theta(i,j) = shrinkageDoubling(&theta_logfcd, X, Lambda, theta,
					     delta0, delta1, Lambda_prior_mean, 
					     Lambda_prior_prec, Lambda_ineq,
					     theta_ineq, *k0, *k1, 
					     *c0, *d0, *c1, 
					     *d1, i, j, z, w, stream, 
					     L, R, param);
	    }
	  }
	}
      }
    
    

      // sample Lambda
      param = 0;
      SampleNoReplace(K, K, K_inds_perm, K_array, stream);   
      for (int ii=0; ii<K; ++ii){
	int i = K_inds_perm[ii];
	SampleNoReplace(D, D, D_inds_perm, D_array, stream);   
	for (int jj=0; jj<D; ++jj){
	  int j = D_inds_perm[jj];
	  if (Lambda_eq(i,j) == -999){
	    w = *lambda_w;
	    p = *lambda_p;
	    L = -1.0;
	    R = 1.0;
	    funval = Lambda_logfcd(Lambda(i,j), X, Lambda, 
				   theta, delta0,
				   delta1, Lambda_prior_mean, 
				   Lambda_prior_prec,
				   Lambda_ineq, theta_ineq, *k0, *k1, *c0, *d0,
				   *c1, *d1, i, j);
	    z = funval - stream.rexp(1.0);
	    if (*method_step == 1){
	      StepOut(&Lambda_logfcd, X, Lambda, theta, delta0, delta1, 
		      Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
		      theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, i, j, z, 
		      w, p, stream, L, R, param);
	      Lambda(i,j) = shrinkageStep(&Lambda_logfcd, X, Lambda, theta,
					  delta0, delta1, Lambda_prior_mean, 
					  Lambda_prior_prec, Lambda_ineq,
					  theta_ineq, *k0, *k1, *c0, *d0, 
					  *c1, *d1, i, j, z, w, stream, 
					  L, R, param);
	    }
	    else{
	      doubling(&Lambda_logfcd, X, Lambda, theta, delta0, delta1, 
		       Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
		       theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, i, j, z, 
		       w, p, stream, L, R, param);
	      Lambda(i,j) = shrinkageDoubling(&Lambda_logfcd, X, Lambda, theta,
					      delta0, delta1, Lambda_prior_mean, 
					      Lambda_prior_prec, Lambda_ineq,
					      theta_ineq, *k0, *k1, *c0, *d0, 
					      *c1, *d1, i, j, z, w, stream, 
					      L, R, param);
	    }
	  }
	}
      }
    

    
      // sample delta0
      param = 2;
      w = *delta0_w;
      p = *delta0_p;
      L = -1.0;
      R =  1.0;    
      funval = delta0_logfcd(delta0, X, Lambda, 
			     theta, delta0,
			     delta1, Lambda_prior_mean, 
			     Lambda_prior_prec,
			     Lambda_ineq, theta_ineq, 
			     *k0, *k1, *c0, *d0,
			     *c1, *d1, 0, 0);
      z = funval - stream.rexp(1.0);
      if (*method_step == 1){
	StepOut(&delta0_logfcd, X, Lambda, theta, delta0, delta1, 
		Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
		theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, 0, 0, z, 
		w, p, stream, L, R, param);
	delta0 = shrinkageStep(&delta0_logfcd, X, Lambda, theta,
			       delta0, delta1, Lambda_prior_mean, 
			       Lambda_prior_prec, Lambda_ineq, theta_ineq,
			       *k0, *k1, *c0, *d0, *c1, *d1, 0, 0,
			       z, w, stream, L, R, param);    
      }
      else{
	doubling(&delta0_logfcd, X, Lambda, theta, delta0, delta1, 
		 Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
		 theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, 0, 0, z, 
		 w, p, stream, L, R, param);
	delta0 = shrinkageDoubling(&delta0_logfcd, X, Lambda, theta,
				   delta0, delta1, Lambda_prior_mean, 
				   Lambda_prior_prec, Lambda_ineq, theta_ineq,
				   *k0, *k1, *c0, *d0, *c1, *d1, 0, 0,
				   z, w, stream, L, R, param);
      }



      // sample delta1
      param = 3;
      w = *delta1_w;
      p = *delta1_p;
      L = -1.0;
      R = 1.0; 
      funval = delta1_logfcd(delta1, X, Lambda, 
			     theta, delta0,
			     delta1, Lambda_prior_mean, 
			     Lambda_prior_prec,
			     Lambda_ineq, theta_ineq, *k0, *k1, *c0, *d0,
			     *c1, *d1, 0, 0);
      z = funval - stream.rexp(1.0);
      if (*method_step == 1){
	StepOut(&delta1_logfcd, X, Lambda, theta, delta0, delta1, 
		Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
		theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, 0, 0, z, 
		w, p, stream, L, R, param);
	delta1 = shrinkageStep(&delta1_logfcd, X, Lambda, theta,
			       delta0, delta1, Lambda_prior_mean, 
			       Lambda_prior_prec, Lambda_ineq, theta_ineq,
			       *k0, *k1, *c0, *d0, *c1, *d1, 0, 0,
			       z, w, stream, L, R, param);
      }
      else{
	doubling(&delta1_logfcd, X, Lambda, theta, delta0, delta1, 
		 Lambda_prior_mean, Lambda_prior_prec, Lambda_ineq,
		 theta_ineq, *k0, *k1, *c0, *d0, *c1, *d1, 0, 0, z, 
		 w, p, stream, L, R, param);
	delta1 = shrinkageDoubling(&delta1_logfcd, X, Lambda, theta,
				   delta0, delta1, Lambda_prior_mean, 
				   Lambda_prior_prec, Lambda_ineq, theta_ineq,
				   *k0, *k1, *c0, *d0, *c1, *d1, 0, 0,
				   z, w, stream, L, R, param);
      }
    
      // print results to screen
      if (verbose[0] > 0 && iter % verbose[0] == 0){
	Rprintf("\n\nMCMCirtKdRob iteration %i of %i \n", (iter+1), tot_iter);
      }
        
      // store results
      if ((iter >= burnin[0]) && ((iter % thin[0]==0))) {      
      
	// store Lambda
	if (storeitem[0]==1){
	  //Matrix<double> Lambda_store_vec = reshape(Lambda,1,K*D);
	  //for (int l=0; l<K*D; ++l)
	  //  Lambda_store(count, l) = Lambda_store_vec[l];
	  rmview(Lambda_store(count, _)) = Lambda;
	}
      
	// store theta
	if (storeability[0]==1){
	  //Matrix<double> theta_store_vec = reshape(theta, 1, N*D);
	  //for (int l=0; l<N*D; ++l)
	  //  theta_store(count, l) = theta_store_vec[l];
	  rmview(theta_store(count, _)) = theta;
	}
      
	// store delta0 and delta1
	delta0_store[count] = delta0;
	delta1_store[count] = delta1;

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



    delete [] K_array;
    delete [] N_array;
    delete [] D_array;
    delete [] Dm1_array;
    delete [] K_inds_perm;
    delete [] N_inds_perm;
    delete [] D_inds_perm;
    delete [] Dm1_inds_perm;



  
    // return output
    Matrix<double> output = delta0_store;
    output = cbind(output, delta1_store);  
    if (*storeitem == 1){
      output = cbind(output, Lambda_store);
    }
    if(*storeability == 1) {
      output = cbind(output, theta_store);
    }

    int size = *samplerow * *samplecol;
    for (int i=0; i<size; ++i)
      sampledata[i] = output[i];
    

}







extern "C"{

  // function called by R to fit model
  void
  irtKdRobpost (double* sampledata, const int* samplerow, 
		const int* samplecol,
		const int* Xdata, const int* Xrow, const int* Xcol,
		const int* burnin, const int* mcmc,  const int* thin,
		const int *uselecuyer, const int *seedarray,
		const int *lecuyerstream, 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* Lamstartdata, const int* Lamstartrow, 
		const int* Lamstartcol, 
		const double* thetstartdata, const int* thetstartrow, 
		const int* thetstartcol, 
		const double* Lameqdata, const int* Lameqrow, 
		const int* Lameqcol,
		const double* Lamineqdata, const int* Lamineqrow, 
		const int* Lamineqcol,
		const double* theteqdata, const int* theteqrow, 
		const int* theteqcol,
		const double* thetineqdata, const int* thetineqrow, 
		const int* thetineqcol,
		const double* Lampmeandata, const int* Lampmeanrow, 
		const int* Lampmeancol,
		const double* Lampprecdata, const int* Lampprecrow,
		const int* Lamppreccol, 
		const double* k0, const double* k1,
		const double* c0, const double* c1,
		const double* d0, const double* d1,
		const int* storeitem,
		const int* storeability
		) {
    
    // put together matrices
    const Matrix<int> X(*Xrow, *Xcol, Xdata);
    Matrix<double> Lambda(*Lamstartrow, *Lamstartcol, Lamstartdata);
    Matrix<double> theta(*thetstartrow, *thetstartcol, thetstartdata);
    const Matrix<double> Lambda_eq(*Lameqrow, *Lameqcol, Lameqdata);
    const Matrix<double> Lambda_ineq(*Lamineqrow, *Lamineqcol, 
				     Lamineqdata);
    const Matrix<double> theta_eq(*theteqrow, *theteqcol, 
				  theteqdata);
    const Matrix<double> theta_ineq(*thetineqrow, *thetineqcol, 
				    thetineqdata);
    const Matrix<double> Lambda_prior_mean(*Lampmeanrow, 
					   *Lampmeancol, 
					   Lampmeandata);
    const Matrix<double> Lambda_prior_prec(*Lampprecrow, 
					   *Lamppreccol, 
					   Lampprecdata);  
    
  
   

    MCMCPACK_PASSRNG2MODEL(MCMCirtKdRob_impl, X, Lambda, theta, 
			   Lambda_eq, Lambda_ineq, theta_eq,
			   theta_ineq, Lambda_prior_mean,
			   Lambda_prior_prec, burnin, mcmc, thin,
			   verbose, 
			   method_step,
			   theta_w, theta_p,
			   lambda_w, lambda_p,
			   delta0_w, delta0_p,
			   delta1_w, delta1_p,
			   delta0start, delta1start,
			   k0,  k1,
			   c0,  c1,
			   d0,  d1,
			   storeitem,
			   storeability,
			   sampledata, samplerow, 
			   samplecol			  
			   );

  
  }

}





#endif