             // 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];
-            // Just assume E(meas|chi) = chi for now
             h_meas[k * n + i] = SED_ik + sigma_msmt[k] * snorm(rng);
 		}
 	}

        double* p_meas = thrust::raw_pointer_cast(&d_meas[0]);
        double* p_meas_unc = thrust::raw_pointer_cast(&d_meas_unc[0]);
        double* p_chi = thrust::raw_pointer_cast(&d_chi[0]);
        double* p_jump_sigma = thrust::raw_pointer_cast(&d_jump_sigma[0]);
        double* p_logdens = thrust::raw_pointer_cast(&d_logdens[0]);
        int current_iter = i + 1;
        update_chi<<<nBlocks,nThreads>>>(p_chi, p_meas, p_meas_unc, n, p_logdens, p_jump_sigma, devStates, current_iter);
        
        // Generate new theta in parallel with GPU calculation above
        double proposed_theta[2];


// Perform the MHA update on the n parameters, done in parallel on the GPU. This is what the kernel does.
__global__
void update_chi(double* chi, double* meas, double* meas_unc, int n, double* logdens, double* prop_cholfact,
                curandState* state, int current_iter)
{
	int i = blockDim.x * blockIdx.x + threadIdx.x;
