             meas_unc_i[k] = meas_unc[n * k + i];
         }
 
-        double logdens_meas = logdensity_meas(meas_i, meas_unc_i, new_chi);
-        double logdens_pop = logdensity_pop(new_chi, c_theta);
-        double logdens_prop = logdens_meas + logdens_pop;
+        double logdens_meas_prop = logdensity_meas(meas_i, meas_unc_i, new_chi);
+        double logdens_pop_prop = logdensity_pop(new_chi, c_theta);
+        double logdens_prop = logdens_meas_prop + logdens_pop_prop;
+
+        bool finite_logdens = isfinite(logdens_prop);
 
         // Compute the Metropolis ratio
-        double ratio = logdens_prop - logdens[i] - logdens_meas;
+        double logdens_old = logdens_pop[i] + logdens_meas[i];
+        double ratio = logdens_prop - logdens_old;
+
+        if (i == 0) {
+			debug_info[0] = logdens_meas_prop - logdens_meas[i];
+			debug_info[1] = logdens_pop_prop - logdens_pop[i];
+			debug_info[2] = ratio;
+		}
+
         ratio = (ratio < 0.0) ? ratio : 0.0;
         ratio = exp(ratio);
 
         // Now randomly decide whether to accept or reject
         double unif_draw = curand_uniform_double(&localState);
 
-        if (unif_draw < ratio) {
+        if ((unif_draw < ratio) && finite_logdens) {
             // Accept this proposal, so save this parameter value and conditional log-posterior
             for (int j=0; j<c_p; j++) {
                 chi[n * j + i] = new_chi[j];
             }
-            logdens[i] = logdens_pop;
+            logdens_pop[i] = logdens_pop_prop;
+            ++naccept[i];
         }
 
         // Copy state back to global memory
         state[i] = localState;
 
-        // Finally, adapt the Cholesky factor of the proposal distribution
-        // TODO: check that the ratio is finite before doing this step
-        double unit_norm = 0.0;
-        for (int j=0; j<c_p; j++) {
-            unit_norm += snorm_deviate[j] * snorm_deviate[j];
-        }
-        unit_norm = sqrt(unit_norm);
-        double decay_sequence = 1.0 / pow((double) current_iter, c_decay_rate);
-        double scaled_coef = sqrt(decay_sequence * fabs(ratio - c_target_rate)) / unit_norm;
-        for (int j=0; j<c_p; j++) {
-            scaled_proposal[j] *= scaled_coef;
-        }
-        bool downdate = (ratio < c_target_rate);
-        CholUpdateR1(cholfactor, scaled_proposal, c_p, downdate);
-
-        // copy cholesky factor for this data point back to global memory
-        cholfact_index = 0;
-        for (int j=0; j<dim_cholfactor; j++) {
-            prop_cholfact[n * j + i] = cholfactor[j];
-        }
-	}
-}
-
+        if (finite_logdens) {
+			// Finally, adapt the Cholesky factor of the proposal distribution
+			double unit_norm = 0.0;
+			for (int j=0; j<c_p; j++) {
+				unit_norm += snorm_deviate[j] * snorm_deviate[j];
+			}
+			unit_norm = sqrt(unit_norm);
+			double decay_sequence = 1.0 / pow((double) current_iter, c_decay_rate);
+			double scaled_coef = sqrt(decay_sequence * fabs(ratio - c_target_rate)) / unit_norm;
+			for (int j=0; j<c_p; j++) {
+				scaled_proposal[j] *= scaled_coef;
+			}
+			bool downdate = (ratio < c_target_rate);
+			CholUpdateR1(cholfactor, scaled_proposal, c_p, downdate);
+
+			// copy cholesky factor for this data point back to global memory
+			for (int j=0; j<dim_cholfactor; j++) {
+				prop_cholfact[n * j + i] = cholfactor[j];
+			}
+        }
+	}
+}
+
+// calculate the logdensity of chi|meas,theta for each chi on the device, needed for computing the
+// log-densities of the initial values
+__global__
+void g_logdens_meas(double* chi, double* meas, double* meas_unc, int ndata, double* logdens_meas)
+{
+	int i = blockDim.x * blockIdx.x + threadIdx.x;
+	if (i < ndata)
+	{
+		double this_chi[p];
+		double this_meas[mfeat];
+		double this_meas_unc[mfeat];
+		for (int j=0; j<p; j++) {
+			// grab the characteristics for this data point
+			this_chi[j] = chi[j * ndata + i];
+		}
+        for (int k=0; k<c_m; k++) {
+            // grab the measurements for this data point
+            this_meas[k] = meas[ndata * k + i];
+            this_meas_unc[k] = meas_unc[ndata * k + i];
+        }
+        logdens_meas[i] = logdensity_meas(this_meas, this_meas_unc, this_chi);
+	}
+}
+
+// calculate the logdensity of theta for each chi on the device, needed for updating theta
+__global__
+void g_logdens_pop(double* chi, int ndata, double* logdens_pop)
+{
+	int i = blockDim.x * blockIdx.x + threadIdx.x;
+	if (i < ndata)
+	{
+		double this_chi[p];
+		for (int j=0; j<p; j++) {
+			this_chi[j] = chi[j * ndata + i];
+		}
+		logdens_pop[i] = logdensity_pop(this_chi, c_theta);
+	}
+}
+
+/*
 struct zsqr : public thrust::unary_function<double,double> {
     double mu, var;
     zsqr(double m, double v) : mu(m), var(v) {}
 
     __device__ __host__
-    double operator()(double chi) {
+    double operator()(double* chi) {
         double chi_cent = chi - mu;
         double logdens_pop = -0.5 * log(var) - 0.5 * chi_cent * chi_cent / var;
         return logdens_pop;
     }
 };
+*/
 
 // Generate a data set and fill the chi, meas, and meas_sigma host vector with the values
 void generate_data(int ndata, thrust::host_vector<double>& theta, thrust::host_vector<double>& chi, thrust::host_vector<double>& meas,

            thetafile << h_theta[j] << " ";
        }
        thetafile << std::endl;
    }
    std::cout << "Number of accepted thetas: " << naccept_theta << std::endl;
    thetafile.close();

    cudaFree(devStates);  // Free up the memory on the GPU from the RNG states

	return 0;
}


        if (unif < ratio) {
            // Accept the proposed theta
            h_theta = proposed_theta;
            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)));
		}

        // finally, adapt the proposal covariance of theta by updating its cholesky factor
        double unit_norm = 0.0;
