                 scaled_proposal_j += cholfactor[cholfact_index] * snorm_deviate[k];
                 cholfact_index++;
             }
-            scaled_proposal_j = 0.0;
+
             new_chi[j] = chi[n * j + i] + scaled_proposal_j;
             scaled_proposal[j] = scaled_proposal_j;
         }


        CUDA_CALL(cudaDeviceSynchronize());
        h_chi = d_chi;
        h_logdens = d_logdens;

        // Compute Metropolis ratio
        //
        // right now we loop over the elements of chi because we assume that they are statistically independent, i.e.,
        // the theta array does not contain correlations among the element of chi. this is unrealistic, and we will need
        // to come up with a better way to do the transform + reduction: maybe use CUBLAS?
        thrust::device_vector<double>::iterator chi_iter_begin = d_chi.begin();

        double logdens_pop = 0.0;

                                                zsqrj, logdens_pop, thrust::plus<double>());
            thrust::advance(chi_iter_begin, ndata);
        }
        */

        // get log-density of current value of theta
        double logdens_old = thrust::reduce(d_logdens_pop.begin(), d_logdens_pop.end());

        // get value of population log-density for proposed theta value for each chi value. note that
        // this overwrites the value of d_logdens_pop.
        double* p_theta = thrust::raw_pointer_cast(&proposed_theta[0]);
        // copy the proposed theta to constant memory on the device
        CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, proposed_theta.size() * sizeof(*p_theta)));
        p_chi = thrust::raw_pointer_cast(&d_chi[0]);  // grab the pointers in case thrust changed the location
        p_logdens_pop = thrust::raw_pointer_cast(&d_logdens_pop[0]);

        g_logdens_pop<<<nBlocks,nThreads>>>(p_chi, ndata, p_logdens_pop);
        CUDA_CALL(cudaDeviceSynchronize());

        double logdens_prop = thrust::reduce(d_logdens_pop.begin(), d_logdens_pop.end());

        double lograt = logdens_prop - logdens_old;
        lograt = std::min(lograt, 0.0);
        double ratio = exp(lograt);
        double unif = uniform(rng);
