Missing data handling

Modeling
11 
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);
}