 // Function to compute the conditional log-density of the measurements given the chi value
 __device__
 double logdensity_meas(double* meas, double* meas_unc, int m, double* chi) {
-    
-    
-    
-    
-    double centered_meas = meas_i - new_chi;
-    double logdens_meas = -0.5 * centered_meas * centered_meas / (meas_unc_i * meas_unc_i);
-
+
+    double normalization, beta, temperature;
+    normalization = exp(chi[0]);
+    beta = chi[1];
+    temperature = exp(chi[2]);
+    
+    double logdensity = 0.0
+    for (int k=0; k<m; k++) {
+        model_sed = modified_bbody(c_nu[k], normalization, beta, temperature);
+        double stnd_meas = (meas[k] - model_sed) / meas_unc[k]
+        logdensity += -0.5 * stnd_meas * stnd_meas;
+    }
+    return logdensity;
 }
 
 // Function to compute the conditional log-density of the chi values given theta
 __device__ __host__
 double logdensity_pop(double* chi, double* theta) {
-    
+    // right now this is just an independent p-dimensional gaussian distribution, for simplicity
+    double logdensity = 0.0;
+    for (int j=0; j<c_p; j++) {
+        double mu_j = theta[j];
+        double var_j = theta[j + c_p];
+        double centered_chi_j = chi[j] - mu_j;
+        logdensity += -0.5 * log(var_j) - 0.5 * centered_chi_j * centered_chi_j / var_j;
+    }
+    return logdensity;
 }
 
 // Initialize the random number generator state

static __constant__ double c_theta[dim_theta];

// Constants used for the model SED
static __constant__ double c_nu0 = 2.3e11;  // nu0 = 230 GHz
static const double nu0 = 2.3e11;
static const double hplanck = 6.6260755e-27;  // Planck's constant, in cgs
static const double clight = 2.997925e10;
static const double kboltz = 1.380658e-16;

	// Create simulated characteristics and data
    thrust::host_vector<double> h_true_chi(n * p);
    double sigma_msmt[5] = {2.2e-4, 3.3e-4, 5.2e-4, 3.7e-4, 2.2e-4};  // Standard deviation for the measurement errors
    
	for (int i=0; i<n; i++) {
        // Loop over the data indices
