 		double* p_sigmas = thrust::raw_pointer_cast(&d_sigmas[0]);
 
 		// cuda grid launch
-		dim3 nThreads(32,8);
+		dim3 nThreads(512,1);
 		dim3 nBlocks((n + nThreads.x-1) / nThreads.x, (n_theta + nThreads.y-1) / nThreads.y);
 		printf("nBlocks: %d  %d\n", nBlocks.x, nBlocks.y);  // no more than 64k blocks!
 		marginals<<<nBlocks,nThreads>>>(p_theta, dim_theta, n_theta, 

	// 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 dmu = (mu_hi - mu_lo)/(n_theta-1.);
	double mu;
	thrust::host_vector<double> h_theta(dim_theta*n_theta);

//#include <cstdlib>

#include <vector>
// constants for integration
#define N_GH 10

