for(j in 1:N) { // see http://www.stata.com/manuals14/rheckman.pdf#rheckmanMethodsandformulas
if (R[j] == 1) {
real err = (y[j] - x[j] * beta) / sigma;
target += normal_lcdf((z[j] * gamma + err * rho) / sqrt(1 - square(rho)) | 0, 1)
- 0.5 * square(err) - log(sqrt(2 * pi()) * sigma);
}
else target += -normal_lcdf(-z[j] * gamma | 0, 1);
}