 	if (i<n)
 	{
 		double mu = theta[0];
-		marg[i] = mu * c_wts[3] * rho(mu);
+		double var = theta[1];
+		double x;
+		double rho;
+		marg[i] = 0;
+		for (int j=0; j<N_GH; j++) {
+			x = features[0][i] + c_rt2*sigmas[0][i]*c_absc[j]
+			rho = exp(-(x-mu)**2/(2.*var))/c_rt2pi
+			marg[i] += c_wts[j]*rho;
+		}
 	}
 }
 	  



// constants for integration
static __constant__ double c_absc[] = {
	-3.4361591188377352, -2.5327316742327906, -1.756683649299881, -1.0366108297895147,
	-0.34290132722370503, 0.34290132722370431, 1.0366108297895129, 1.7566836492998816,

{
	// measurements
	int n = 10; // # of items
	int m = 1; // # of features
	wrapvec d_features(m,n);
	wrapvec d_sigmas(m,n);

	unsigned int seed = 9872498732;
	static thrust::minstd_rand rng(seed);
	thrust::random::experimental::normal_distribution<double> dist(0.0, 1.0);
	// thrust::generate(d_vec.begin(), d_vec.end(), dist(rng));

	for (int i=0; i<n; i++) {
		for (int j=0; j<m; j++) {
			d_features[j][i] = 5.;
			d_sigmas[j][i] = 1.;
		}
	}

	// alloc mem for marginals
	thrust::device_vector<double> d_marg(n);
