     // Initialize the chi values, as well as the cholesky factors for the theta and chi proposals
     initialize_parameters(ndata, h_theta, h_chi, h_theta_cholfact, h_chi_cholfactor, h_true_chi);
 
-
     // Allocate memory for arrays on device and copy the values from the host
     thrust::device_vector<double> d_meas = h_meas;
     thrust::device_vector<double> d_meas_unc = h_meas_unc;
     thrust::device_vector<double> d_chi = h_chi;
     thrust::device_vector<double> d_chi_cholfactor = h_chi_cholfactor;
-    thrust::device_vector<double> d_logdens(ndata);  // log posteriors for an individual data point
-
-    // Initialize the conditional log-posterior of chi|theta
-    thrust::fill(d_logdens.begin(), d_logdens.end(), -1.0e300);
-
-	// Load the initial guess for theta into constant memory
+    thrust::device_vector<double> d_logdens_meas(ndata);  // log densities of chi values
+    thrust::device_vector<double> d_logdens_pop(ndata); // log densities of theta values for each chi
+
+   	// Load the initial guess for theta into constant memory
     double* p_theta = &h_theta[0];
 	CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, h_theta.size() * sizeof(*p_theta)));
 


// Global random number generator and distributions for generating random numbers on the host. The random number generator used
// is the Mersenne Twister mt19937 from the BOOST library.

boost::random::mt19937 rng;
boost::random::normal_distribution<> snorm(0.0, 1.0); // Standard normal distribution
boost::random::uniform_real_distribution<> uniform(0.0, 1.0); // Uniform distribution from 0.0 to 1.0

void initialize_parameters(int ndata, thrust::host_vector<double>& theta, thrust::host_vector<double>& chi,
		thrust::host_vector<double>& theta_cholfactor, thrust::host_vector<double>& chi_cholfactor, thrust::host_vector<double>& true_chi) {

	// For now, initialize chi values for drawing from their parent distribution.
	for (int i = 0; i < ndata; ++i) {
		for (int j = 0; j < p; ++j) {
			double chi_mean = theta[j];
			double chi_sig = sqrt(exp(theta[p + j]));
			chi[j * ndata + i] = chi_mean + chi_sig * snorm(rng);
		}
	}

	for (int j = 0; j < dim_theta; ++j) {
		theta[j] *= 0.2 * snorm(rng);
	}

	// Initialize the cholesky factor for the theta proposals
    int diag_index = 0;
    for (int j=0; j<dim_theta; j++) {
    	// Initialize the covariance matrix of the theta proposal distribution to be 0.01 * identity(p)
        theta_cholfactor[diag_index] = 0.01;
        diag_index += j + 2;
    }

