         for (int j=0; j<dim_theta; j++) {
             thetafile << h_theta[j] << " ";
         }
-        thetafile << std::endl;
+        double logdens_pop = thrust::reduce(d_logdens_pop.begin(), d_logdens_pop.end());
+        thetafile << logdens_pop << std::endl;
         // Save the chi values
         for (int j=0; j<p; j++) {
         	chifile << h_chi[ndata * j] << " ";
         }
-        chifile << std::endl;
+        double logdens_meas = thrust::reduce(d_logdens_meas.begin(), d_logdens_meas.end());
+        chifile << logdens_meas << std::endl;
     }
     thetafile.close();
 

        return 2;
    }

    /*
    Set up parallel random number generators on the GPU
    */
    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 *
                         sizeof(curandState)));
    // Initialize the random number generator states on the GPU
    setup_kernel<<<nBlocks,nThreads>>>(devStates);

    // Wait until RNG stuff is done running on the GPU, make sure everything went OK
    CUDA_CALL(cudaDeviceSynchronize());

    /*******************************
    	Now run the MCMC sampler
	*******************************/

    //std::ofstream chifile("chis.dat");
    std::ofstream thetafile("thetas.dat");
    std::ofstream chifile("chi0.dat");
    int mcmc_iter = 25000;
    int naccept_theta = 0;
    std::cout << "Running MCMC Sampler...." << std::endl;


                chi[n * j + i] = new_chi[j];
            }
            logdens_pop[i] = logdens_pop_prop;
            logdens_meas[i] = logdens_meas_prop;
            ++naccept[i];
        }

