     double logdensity = 0.0;
     for (int j=0; j<dim_chi; j++) {
         double mu_j = theta[j];
-        double var_j = exp(theta[j + ntheta]);
+        double var_j = exp(theta[j + dim_chi]);
         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;
+}
+
+// Function to compute the conditional log-density of the chi values given theta
+double h_logdensity_pop(double* chi, double* theta) {
+    int dim_chi = p;
+    int ntheta = dim_theta;
+    // right now this is just an independent p-dimensional gaussian distribution, for simplicity
+    std::cout << "theta: ";
+    for (int j=0; j<dim_theta; j++) {
+        std::cout << " " << theta[j];
+    }
+    std::cout << std::endl;
+    double logdensity = 0.0;
+    for (int j=0; j<dim_chi; j++) {
+        double mu_j = theta[j];
+        double var_j = exp(theta[j + dim_chi]);
+        double centered_chi_j = chi[j] - mu_j;
+        std::cout << "mu_j: " << mu_j << ", var_j: " << var_j << ", chi[j]: " << chi[j] << std::endl;
         logdensity += -0.5 * log(var_j) - 0.5 * centered_chi_j * centered_chi_j / var_j;
     }
     return logdensity;

	}
}

struct zsqr {
    double mu, var;
    zsqr(double m, double v) : mu(m), var(v) {}
    
    __device__ __host__
    double operator()(double chi) {
        double chi_cent = chi - mu;
        double chisqr = -0.5 * chi_cent * chi_cent / var;
        return chisqr;
    }
};


        }
        
        CUDA_CALL(cudaDeviceSynchronize());
        h_chi = d_chi;
        h_logdens = d_logdens;
        
        h_mha_ratios = d_mha_ratios;
        double min_mha_ratio = 1e10;
