     h_theta[4] = log(sig_beta * sig_beta);
     h_theta[5] = log(sig_temp * sig_temp);
 
-    // Construct initial cholesky factor for the covariance matrix of the theta proposal distribution
+    // Allocate memory for arrays of chi, measurement, and measurement uncertainties on host
+    thrust::host_vector<double> h_meas(ndata * mfeat);  // The measurements, mfeat values for each of ndata data points
+	thrust::host_vector<double> h_meas_unc(ndata * mfeat);  // The measurement uncertainties
+    thrust::host_vector<double> h_chi(ndata * p);  // Unknown characteristics, p values for each of ndata data points
+    thrust::host_vector<double> h_true_chi(ndata * p);
+
+    // Generate characteristics and data
+    thrust::host_vector<double> h_true_chi(ndata * p); // true values of the chi-values
+    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
+
+    std::cout << "Generating some data..." << std::endl;
+    generate_data(ndata, h_theta, h_true_chi, h_meas, h_meas_unc);
+
+    // Allocate memory for cholesky factors on host.
+    //
+    // Cholesky factor of chi proposal distribution for each data point, used by the Metropolis algorithm. The proposal
+    // covariance matrix of chi for each data point is Covar_i = PropChol_i * transpose(PropChol_i).
+    thrust::host_vector<double> h_chi_cholfactor(ndata * dim_cholfactor);
+    thrust::fill(h_chi_cholfactor.begin(), h_chi_cholfactor.end(), 0.0);
+
+    // Cholesky factor for the covariance matrix of the theta proposal distribution
     int dim_cholfactor_theta = dim_theta * dim_theta - ((dim_theta - 1) * dim_theta) / 2;
     thrust::host_vector<double> h_theta_cholfact(dim_cholfactor_theta);
-    int diag_index = 0;
-    for (int j=0; j<dim_theta; j++) {
-        h_theta_cholfact[diag_index] = 1e-2;
-        diag_index += j + 2;
-    }
-
-    // Allocate memory for arrays on host
-    thrust::host_vector<double> h_meas(ndata * mfeat);  // The measurements, m values for each of n data points
-	thrust::host_vector<double> h_meas_unc(ndata * mfeat);  // The measurement uncertainties
-    thrust::host_vector<double> h_chi(ndata * p);  // Unknown characteristics, p values for each of n data points
-
-    // 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).
-    thrust::host_vector<double> h_prop_cholfact(ndata * dim_cholfactor);
-    thrust::fill(h_prop_cholfact.begin(), h_prop_cholfact.end(), 0.0);
-
-	// Create simulated characteristics and data
-    thrust::host_vector<double> h_true_chi(ndata * 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
-
-    std::cout << "Generating some data..." << std::endl;
-
-    // Initialize the chis to their true values for now
-    thrust::copy(h_true_chi.begin(), h_true_chi.end(), h_chi.begin());
+
+    // Initialize the chi values, as well as the cholesky factors for the theta and chi proposals
+    initialize_parameters(ndata, h_theta, h_chi, h_theta_cholfact, h_chi_cholfactor, h_true_chi);
+
 
     // Allocate memory for arrays on device and copy the values from the host
     thrust::device_vector<double> d_meas = h_meas;
     thrust::device_vector<double> d_meas_unc = h_meas_unc;
     thrust::device_vector<double> d_chi = h_chi;
-    thrust::device_vector<double> d_prop_cholfact = h_prop_cholfact;
+    thrust::device_vector<double> d_chi_cholfactor = h_chi_cholfactor;
     thrust::device_vector<double> d_logdens(ndata);  // log posteriors for an individual data point
 
+    // Initialize the conditional log-posterior of chi|theta
     thrust::fill(d_logdens.begin(), d_logdens.end(), -1.0e300);
 
-	// Load a single set of thetas into constant memory
+	// Load the initial guess for theta into constant memory
     double* p_theta = &h_theta[0];
 	CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, h_theta.size() * sizeof(*p_theta)));
 

        return 2;
    }

    curandState* devStates;  // Create state object for random number generators on the GPU
    // Allocate memory on GPU for RNG states
    CUDA_CALL(cudaMalloc((void **)&devStates, nThreads.x * nBlocks.x *


	for (int i=0; i<ndata; i++) {
        // Loop over the data indices
		for (int j=0; j<p; j++) {
            // Loop over the chi indices to generate true value of the characteristics
			chi[i + ndata * j] = theta[j] + sqrt(exp(theta[p + j])) * snorm(rng);
        }
        for (int k=0; k<mfeat; k++) {
            // Loop over the feature indices to generate measurements, given the characteristics
