functions { vector sqrt_vec(vector x) { vector[dims(x)[1]] res; for (m in 1:dims(x)[1]){ res[m] = sqrt(x[m]); } return res; } } data { int P; int N; matrix[N,P] X; vector[N] threshvar; vector[N] y; real minthresh; real maxthresh; real vary; } parameters { real muL; real muH; vector[P] arLpre; vector[P] arHpre; real thresh; real s2L; real s2H; real g1a; real g1b; real g2a; real g2b; vector[P] l1a; vector[P] l1b; vector[P] l2a; vector[P] l2b; } transformed parameters { vector[N] trans; vector[P] arL; vector[P] arH; real g1; real g2; vector[P] l1; vector[P] l2; g1=g1a*sqrt(g1b); g2=g2a*sqrt(g2b); l1=l1a .* sqrt_vec(l1b); l2=l2a .* sqrt_vec(l2b); arL=arLpre .* (l1*g1); arH=arHpre .* (l2*g2); for(k in 1:N){ if(threshvar[k]