         if (unif < ratio) {
             // Accept the proposed theta
             h_theta = proposed_theta;
-            // grab the pointer in case thrust changed the memory location on us
-            double* p_theta = thrust::raw_pointer_cast(&h_theta[0]);
-            // copy the new theta to constant memory on the device
+            h_logdens_pop = d_logdens_pop;
+            naccept_theta++;
+        } else {
+			// keep current value of theta, but need to restore the values of d_logdens_pop overwritten by
+        	// g_logdens_pop and move the current value of theta back to constant memory on the device
+        	d_logdens_pop = h_logdens_pop;
+            p_theta = thrust::raw_pointer_cast(&h_theta[0]);
+            // copy the proposed theta to constant memory on the device
             CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, h_theta.size() * sizeof(*p_theta)));
-            naccept_theta++;
-        }
+		}
 
         // finally, adapt the proposal covariance of theta by updating its cholesky factor
         double unit_norm = 0.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) {

	// Initalize the chi values to their true values for now
	thrust::copy(true_chi.begin(), true_chi.end(), chi.begin());

	// 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] = 1e-2;
        diag_index += j + 2;
    }



// Perform the MHA update on the n characteristics, done in parallel on the GPU.
__global__
void update_chi(double* chi, double* meas, double* meas_unc, int n, double* logdens_meas, double* logdens_pop,
		double* prop_cholfact, curandState* state, int current_iter, int* naccept, double* debug_info)
{
	int i = blockDim.x * blockIdx.x + threadIdx.x;
	if (i < n)
