             for (int k=0; k<(j+1); k++) {
                 // transform the unit proposal to the centered proposal, drawn from a multivariate normal
                 // or t-distribution.
-
-                // DOUBLE CHECK THIS: DOES NOT LOOK RIGHT, WHY DOES K NOT APPEAR????
-                scaled_proposal_j += cholfactor[cholfact_index] * snorm_deviate[j];
+                scaled_proposal_j += cholfactor[cholfact_index] * snorm_deviate[k];
                 cholfact_index++;
             }
             scaled_proposal_j = 0.0;

        h_chi = d_chi;
        h_logdens = d_logdens;

        h_mha_ratios = d_mha_ratios;
        double min_mha_ratio = 1e10;
        double max_mha_ratio = -1.0;
        for (int idata=0; idata<n; idata++) {
            min_mha_ratio = std::min(min_mha_ratio, h_mha_ratios[idata]);
            max_mha_ratio = std::max(max_mha_ratio, h_mha_ratios[idata]);
        }
        std::cout << std::endl;
        std::cout << "Iteration: " << current_iter << std::endl;
        std::cout << "Minimum, Maximum MHA ratio: " << min_mha_ratio << ", " << max_mha_ratio << std::endl;

        /*
        std::cout << "Proposed theta: ";
        for (int j=0; j<dim_theta; ++j) {
            std::cout << proposed_theta[j] << " ";
        }
        std::cout << std::endl;
        std::cout << "Current theta: ";
        for (int j=0; j<dim_theta; ++j) {
            std::cout << h_theta[j] << " ";
        }
        std::cout << std::endl;
         */
        // 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?

        double* p_chi = thrust::raw_pointer_cast(&d_chi[0]);
        double* p_prop_cholfact = thrust::raw_pointer_cast(&d_prop_cholfact[0]);
        double* p_logdens = thrust::raw_pointer_cast(&d_logdens[0]);
        int current_iter = i + 1;
        update_chi<<<nBlocks,nThreads>>>(p_chi, p_meas, p_meas_unc, n, p_logdens, p_prop_cholfact, devStates, current_iter);

        // Generate new theta in parallel with GPU calculation above
        thrust::host_vector<double> proposed_theta(dim_theta);
