     // 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)));
 


        CUDA_CALL(cudaDeviceSynchronize());
        h_chi = d_chi;
        h_logdens = d_logdens;

        // Compute Metropolis ratio
        //
        // right now we loop over the elements of chi because we assume that they are statistically independent, i.e.,
        // the theta array does not contain correlations among the element of chi. this is unrealistic, and we will need
        // to come up with a better way to do the transform + reduction: maybe use CUBLAS?
        thrust::device_vector<double>::iterator chi_iter_begin = d_chi.begin();

        double logdens_pop = 0.0;


int main(void)
{
    /* TODO: Add command line options */

    /*
     Population level parameters: theta, where theta parameterizes the distribution of chi.
