 		double x;
 		double r;
 		double mrg = 0;
-		double pi = features[i]; // here only 1 feature!!!
-		double si = sigmas[i];
+		double mi = meas[i]; // here only 1 feature!!!
+		double ui = meas_unc[i];
 		for (int j=0; j<N_GH; j++) {
-			x = pi + c_rt2 * si * c_absc[j];
+			x = mi + c_rt2 * ui * c_absc[j];
 			r = rho(x, mu, var);
 			mrg += c_wts[j] * r;
 		}

	// Create simulated data; copies to GPU dumbly!
	for (int i=0; i<n; i++) {
		for (int j=0; j<m; j++) {
			h_features[i*m+0] = mu_popn_true + dist(rng);
			h_sigmas[i*m+1] = sigma_popn_true;
		}
	}

	// Copy data to GPU
	thrust::device_vector<double> d_features = h_features;
	thrust::device_vector<double> d_sigmas = h_sigmas;

	// Create array of hyperparameter values.
	// Here we condition on sigma_popn_true (for simple analytical result).

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

	// Allocate data on CPU and load 
	thrust::host_vector<double> h_meas(n*m);
	thrust::host_vector<double> h_meas_unc(n*m);

	unsigned int seed = 98724732;
	double sigma_msmt = 3.;
