 	// Currently ineffecient; should build host vector and copy over.
 	int dim_theta = 2;
 	int n_theta = 11;
-	double mu_lo = mu_popn_true - 2*sigma_msmt/sqrt(n);
-	double mu_hi = mu_popn_true + 2*sigma_msmt/sqrt(n);
+	double sqrt_n = sqrt((double)n);
+	double mu_lo = mu_popn_true - 2*sigma_msmt/sqrt_n;
+	double mu_hi = mu_popn_true + 2*sigma_msmt/sqrt_n;
 	double dmu = (mu_hi - mu_lo)/(n_theta-1.);
 	double mu;
 	thrust::host_vector<double> h_theta(dim_theta*n_theta);

//#include <cstdlib>

#include <vector>
/*
typedef thrust::device_vector<double> dvec;
typedef std::vector<dvec> vdvec;

struct wrapvec
{
	vdvec v;

	wrapvec(int m, int n) : v(m)
	{
		for (int i=0; i<m; i++)	{			
			v[i].reserve(n);
		}
	}

	double** ptrs()
	{
		thrust::host_vector<double*> h_ptr(v.size());
		for (unsigned int i=0; i<v.size(); i++)	
			h_ptr[i] = (double*) thrust::raw_pointer_cast(&(v[i][0]));

		thrust::device_vector<double*> d_ptr = h_ptr;
		return (double**) thrust::raw_pointer_cast(&d_ptr[0]); 
	}
};
*/
// constants for integration
#define N_GH 10


}
	  

int main(int argc, char *argv[])
{
	// measurements
	int n = 10000000; // # of items
	int m = 1; // # of features

	int devId = 0;

	if (argc > 1) n = atoi(argv[1]);
	if (argc > 2) devId = atoi(argv[2]);

	cudaError_t err = cudaSetDevice(devId);
	if (err != CUDA_SUCCESS) { return 1; }


	thrust::host_vector<double> h_features(n*m);
	thrust::host_vector<double> h_sigmas(n*m);

