         double mu_j = theta[j];
         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;

        int cholfact_index = 0;
        for (int j=0; j<dim_theta; j++) {
            for (int k=0; k<(j+1); k++) {
                // CHECK THIS: SHOULDN'T THIS BE SNORM_DEVIATE[K]?????
                scaled_proposal[j] += h_theta_cholfact[cholfact_index] * snorm_deviate[j];
                cholfact_index++;
            }
            scaled_proposal[j] = 0.0;

        for (int k=0; k<m; k++) {
            // Loop over the feature indices to generate measurements, given the characteristics
            double bbody_numer = 2.0 * hplanck * nu[k] * nu[k] * nu[k] / (clight * clight);
            double temperature = exp(h_true_chi[n * 2 + i]); // Grab the temperature for this data point
            double bbody_denom = exp(hplanck * nu[k] / (kboltz * temperature)) - 1.0;
            double bbody = bbody_numer / bbody_denom;
            // Grab the SED normalization and power-law index
