     
     // Cholesky factor of proposal distribution for each data point, used by the Metropolis algorithm. The proposal
     // covariance matrix for each data point is Covar_i = PropChol_i * transpose(PropChol_i).
-    dim_cholfactor = p * p - ((p - 1) * p) / 2; // cholesky factor is a lower-diagonal matrix
     thrust::host_vector<double> h_prop_cholfact(n * dim_cholfactor);
     thrust::fill(h_prop_cholfact.begin(), h_prop_cholfact.end(), 0.0);
     

    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
        int diag_index = 0
		for (int j=0; j<p; j++) {
            // Loop over the chi indices to generate true value of the characteristics
			h_true_chi[i + n * j] = h_theta[j] + sqrt(exp(h_theta[p + j])) * snorm(rng);
            // Initialize the covariance matrix of the chi proposal distribution to be 0.01 * identity(p)
            h_prop_cholfact[i + n * diag_index[j]] = 0.01;
            diag_index += j + 2;
        }
        for (int k=0; k<m; k++) {
            // Loop over the feature indices to generate measurements, given the characteristics
            double bbody_numer = 2.0 * hplanck * frequencies[k] * frequencies[k] * frequencies[k] / (clight * clight);
            double temperature = exp(h_true_chi[n * 2 + i]); // Grap the temperature for this data point
            double bbody_denom = exp(hplanck * frequencies[k] / (kboltz * temperature)) - 1.0;
            double bbody = bbody_numer / bbody_denom;
            // Grab the SED normalization and power-law index
            double normalization = exp(h_true_chi[i]);
            double beta = h_true_chi[n + i];
            // Compute the model SED
            double nu_over_nu0 = frequencies[k] / nu0
            double SED_ik = normalization * pow(nu_over_nu0, beta) * bbody;
            // Generate measured SEDs
			h_meas_unc[k * n + i] = sigma_msmt[k];

        }

        // Compute the conditional log-posterior of the proposed parameter values for this data point
        double meas_i[m];
        double meas_unc_i[m];
        for (int k=0; k<c_m; k++) {
            // grab the measurements for this data point
            meas_i[k] = meas[n * k + i];
            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;
        
        // Compute the Metropolis ratio
