 		double* p_logdens_meas = thrust::raw_pointer_cast(&d_logdens_meas[0]);
 		double* p_logdens_pop = thrust::raw_pointer_cast(&d_logdens_pop[0]);
 
+		// Allocate memory on GPU for RNG states
+		CUDA_CHECK_RETURN(cudaMalloc((void **)&p_devStates, nThreads.x * nBlocks.x * sizeof(curandState)));
+		// Initialize the random number generator states on the GPU
+		InitializeRNG<<<nBlocks,nThreads>>>(p_devStates);
+
+		// Wait until RNG stuff is done running on the GPU, make sure everything went OK
+		CUDA_CHECK_RETURN(cudaDeviceSynchronize());
+
 		// set initial values for the characteristics. this will launch a CUDA kernel.
 		InitialValue<ChiType><<<nBlocks,nThreads>>>(p_chi, p_meas, p_meas_unc, p_cholfact, p_logdens_meas,
 				p_logdens_pop, ndata, mfeat, pchi);

	int pchi = 2;

public:
	// Constructor when storing the measurements in std::vector
	DataAugmentation(vecvec& meas, vecvec& meas_unc, dim3& nB, dim3& nT) : nBlocks(nB), nThreads(nT)
	{
		int size1 = meas.size();
		int size2 = meas[0].size();

		ndata = max(size1,size2); // assume that ndata < mfeat and figure out the values from the size
		mfeat = min(size1,size2); // of the arrays.

		_SetArraySizes();

		// copy input data to data members
		for (int j = 0; j < mfeat; ++j) {
			for (int i = 0; i < ndata; ++i) {
				if (size1 < size2) {
					h_meas[ndata * j + i] = meas[j][i];
					h_meas_unc[ndata * j + i] = meas_unc[j][i];
				} else {
					h_meas[ndata * j + i] = meas[i][j];
					h_meas_unc[ndata * j + i] = meas_unc[i][j];
				}
			}
		}
		// copy data from host to device
		d_meas = h_meas;
		d_meas_unc = h_meas_unc;

		thrust::fill(h_cholfact.begin(), h_cholfact.end(), 0.0);
		d_cholfact = h_cholfact;

	// cholesky factors of Metropolis proposal covariance matrix
	hvector h_cholfact;
	dvector d_cholfact;
	// state of parallel random number generator on the device
	curandState* p_devStates;
	// CUDA kernel launch specifications
	dim3& nBlocks;
	dim3& nThreads;
	// MCMC sampler parameters
	int current_iter;
	double target_rate;
	double decay_rate;
	thrust::host_vector<int> h_naccept;
	thrust::device_vector<int> d_naccept;
};

// Base class for a population level parameter
