     //std::ofstream chifile("chis.dat");
     std::ofstream thetafile("thetas.dat");
     
-    int mcmc_iter = 10;
+    int mcmc_iter = 2;
     int naccept_theta = 0;
     std::cout << "Running MCMC Sampler...." << std::endl;
     
     thrust::host_vector<double> h_mha_ratios(n);
     thrust::device_vector<double> d_mha_ratios(n);
+    thrust::host_vector<double> h_logdens;
     
     for (int i=0; i<mcmc_iter; i++) {
         // Now grab the pointers to the vectors, needed to run the kernel since it doesn't understand Thrust

            min_mha_ratio = std::min(min_mha_ratio, h_mha_ratios[idata]);
            max_mha_ratio = std::max(max_mha_ratio, h_mha_ratios[idata]);
        }
        std::cout << "Iteration: " << current_iter << std::endl;
        std::cout << "Minimum, Maximum MHA ratio: " << min_mha_ratio << ", " << max_mha_ratio << std::endl;
        

    thrust::fill(d_logdens.begin(), d_logdens.end(), -1.0e300);
    
	// Load a single set of thetas into constant memory
    double* p_theta = &h_theta[0];
	CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, h_theta.size() * sizeof(*p_theta)));
    
    // Cuda grid launch
