Dear stan users:
I write a likelihood function and use optimizing function to find MLE for the parameters (which are bivariate normal mean, rho and sigma). I wanted to test on different combinations of N ( the number of classical datapoints) , Rho ( the correlation of the underlying data) and M (the number of “symbols” constructed from these underlying classical data and my model only uses these symbols to estimate the bivariate normal distribution parameters rather than using the classical data points).
I have tested the model, it works fine. However, when I want to test on the following 45 combinations (N, RHO, M). The R session freezes for few seconds and pops up r encountered a fatal error the session was terminated . However, when I tried selected combination, say the 45th one, it works fine and MLE looks reasonable, while I tried say 20th combination it experiences the error. so I think my model construction is okay in this case, but it seems that it cannot handle particular cases or…?
I am not sure whether it has something to do with my data combination that the model just cannot handle or it has something to do with rstan. Any advice would be truly appreciated.
The toy example datasets are following
define varying parameters
n=c(60,80,100) ; r = c(0,0.3,0.5,0.7,0.9) ; m = c(20,40,50)
combo = expand.grid(n, r, m)
sigma1 = sigma2= 0.5
set.seed(100)
interval_x = interval_y = index = vector(“list”, length = nrow(combo))
for(i in 1: nrow(combo)){
interval_x[[i]] = interval_y[[i]]= matrix(0, combo$Var3[i], 2)
index[[i]]= matrix(0, combo$Var3[i], 4)
for(j in 1: combo$Var3[i]){
aa = mvrnorm(combo$Var1[i],c(2,5), Sigma = matrix(c(sigma1^2, combo$Var2[i]*sigma1 * sigma2, combo$Var2[i]*sigma1 * sigma2,
sigma2^2),2,2))
index[[i]][j,] = c(which.max(aa[,1]), which.max(aa[,2]),
which.min(aa[,1]), which.min(aa[,2]))
interval_x[[i]][j,]<-c(min(aa[,1]),
max(aa[,1]))
interval_y[[i]][j,]<-c(min(aa[,2]),
max(aa[,2]))
}
}
record how many symbols in each combo are in (2, 21, 3, 31,32, 33 and 4-point case) and record the corresponding intervals
there are altogether 7 cases
cases = vector(“list”, nrow(combo))
interval_x_two = interval_y_two = interval_x_two1 = interval_y_two1=
interval_x_three = interval_y_three =interval_x_three1 = interval_y_three1=interval_x_three2 =
interval_y_three2=interval_x_three3 = interval_y_three3 = interval_x_four = interval_y_four = vector(“list”, nrow(combo))
for(i in 1: nrow(combo)){
cases[[i]]= rep(0,7)
##2-point topright and bottomleft case
cases[[i]][1] = length(which(index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] ))
interval_x_two[[i]] = interval_y_two[[i]] = matrix(0, cases[[i]][1], 2)
interval_x_two[[i]] = interval_x[[i]][which(index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] ),]
interval_y_two[[i]] = interval_y[[i]][which(index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] ),]
##2-point topleft and bottomright case
cases[[i]][2] = length(which(index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] ))
interval_x_two1[[i]] = interval_y_two1[[i]] = matrix(0, cases[[i]][2], 2)
interval_x_two1[[i]] = interval_x[[i]][which(index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] ),]
interval_y_two1[[i]] = interval_y[[i]][which(index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] ),]
##3-point topright case
cases[[i]][3] = length(setdiff(which(index[[i]][,1] == index[[i]][,2]) , which (index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] )))
interval_x_three[[i]] = interval_y_three[[i]] = matrix(0, cases[[i]][3], 2)
interval_x_three[[i]] = interval_x[[i]][setdiff(which(index[[i]][,1] == index[[i]][,2]) , which (index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] )),]
interval_y_three[[i]] = interval_y[[i]][setdiff(which(index[[i]][,1] == index[[i]][,2]) , which (index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] )),]
##3-point bottomleft case
cases[[i]][4] = length(setdiff(which(index[[i]][,3] == index[[i]][,4]) , which (index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] )))
interval_x_three1[[i]] = interval_y_three1[[i]] = matrix(0, cases[[i]][4], 2)
interval_x_three1[[i]] = interval_x[[i]][setdiff(which(index[[i]][,3] == index[[i]][,4]) , which (index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] )),]
interval_y_three1[[i]] = interval_y[[i]][setdiff(which(index[[i]][,3] == index[[i]][,4]) , which (index[[i]][,1] == index[[i]][,2] & index[[i]][,3] == index[[i]][,4] )),]
##3-point bottomright case
cases[[i]][5] = length (setdiff(which(index[[i]][,1] == index[[i]][,4]) ,which(index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] ) ))
interval_x_three2[[i]] = interval_y_three2[[i]] = matrix(0, cases[[i]][5], 2)
interval_x_three2[[i]] = interval_x[[i]][setdiff(which(index[[i]][,1] == index[[i]][,4]) , which (index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] )),]
interval_y_three2[[i]] = interval_y[[i]][setdiff(which(index[[i]][,1] == index[[i]][,4]) , which (index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] )),]
##3-point topleft case
cases[[i]][6] = length(setdiff(which(index[[i]][,2] == index[[i]][,3]) ,which(index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] ) ))
interval_x_three3[[i]] = interval_y_three3[[i]] = matrix(0, cases[[i]][6], 2)
interval_x_three3[[i]] = interval_x[[i]][setdiff(which(index[[i]][,2] == index[[i]][,3]) ,which(index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] ) ),]
interval_y_three3[[i]] = interval_y[[i]][setdiff(which(index[[i]][,2] == index[[i]][,3]) ,which(index[[i]][,2] == index[[i]][,3] & index[[i]][,1] == index[[i]][,4] ) ),]
##4-point case
cases[[i]][7] = combo$Var3[i] - sum(cases[[i]][1:6])
interval_x_four[[i]] = interval_y_four[[i]] = matrix(0, cases[[i]][7], 2)
interval_x_four[[i]] = interval_x[[i]][ which(index[[i]][,1] != index[[i]][,2] & index[[i]][,1] != index[[i]][,3] & index[[i]][,1] != index[[i]][,4] &
index[[i]][,2] != index[[i]][,3] & index[[i]][,2] != index[[i]][,4] &
index[[i]][,3] != index[[i]][,4]),]
interval_y_four[[i]] = interval_y[[i]][ which(index[[i]][,1] != index[[i]][,2] & index[[i]][,1] != index[[i]][,3] & index[[i]][,1] != index[[i]][,4] &
index[[i]][,2] != index[[i]][,3] & index[[i]][,2] != index[[i]][,4] &
index[[i]][,3] != index[[i]][,4]),]
}
The Rstan code as follows
functions{
real binormal_cdf(real z1, real z2, real rho) {
if (z1 != 0 || z2 != 0) {
real denom = fabs(rho) < 1.0 ? sqrt((1 + rho) * (1 - rho)) : not_a_number();
real a1 = (z2 / z1 - rho) / denom;
real a2 = (z1 / z2 - rho) / denom;
real product = z1 * z2;
real delta = product < 0 || (product == 0 && (z1 + z2) < 0);
return 0.5 * (Phi(z1) + Phi(z2) - delta) - owens_t(z1, a1) - owens_t(z2, a2);
}
return 0.25 + asin(rho) / (2 * pi());
}
}
data {
int<lower=0> n;//total number of classical data points per symbol
int<lower=0> m[7];// number of symbols in each of 2,3,4 cases m[1] 2-point case 1
//m[2] 2-point case2; m[3] 3-point bottomleft, m[4] 3-point bottomright, m[5] 3-point topright m[6] 3-point topleft
//m[7] 4-point
vector[2] intervals_x_two[m[1]];//matrix of bivaraite min max in x
vector[2] intervals_y_two[m[1]];//matrix of bivaraite min max in y
vector[2] intervals_x_two1[m[2]];//matrix of bivaraite min max in x
vector[2] intervals_y_two1[m[2]];//matrix of bivaraite min max in y
vector[2] intervals_x_three[m[3]];//matrix of bivaraite min max in x
vector[2] intervals_y_three[m[3]];//matrix of bivaraite min max in y
vector[2] intervals_x_three1[m[4]];//matrix of bivaraite min max in x
vector[2] intervals_y_three1[m[4]];//matrix of bivaraite min max in y
vector[2] intervals_x_three2[m[5]];//matrix of bivaraite min max in x
vector[2] intervals_y_three2[m[5]];//matrix of bivaraite min max in y
vector[2] intervals_x_three3[m[6]];//matrix of bivaraite min max in x
vector[2] intervals_y_three3[m[6]];//matrix of bivaraite min max in y
vector[2] intervals_x[m[7]];//matrix of bivaraite min max in x 4 points
vector[2] intervals_y[m[7]];//matrix of bivaraite min max in y 4 points
}
transformed data{
// for 2 points case 1
vector[2] topright2[m[1]];
vector[2] bottomright2[m[1]];
vector[2] topleft2[m[1]];
vector[2] bottomleft2[m[1]];
//for 2 points case 2
vector[2] topright21[m[2]];
vector[2] bottomright21[m[2]];
vector[2] topleft21[m[2]];
vector[2] bottomleft21[m[2]];
//for three points case 1
vector[2] topright_3[m[3]];
vector[2] bottomright_3[m[3]];
vector[2] topleft_3[m[3]];
vector[2] bottomleft_3[m[3]];
//for three points case 2
vector[2] topright_31[m[4]];
vector[2] bottomright_31[m[4]];
vector[2] topleft_31[m[4]];
vector[2] bottomleft_31[m[4]];
//for three points case 3
vector[2] topright_32[m[5]];
vector[2] bottomright_32[m[5]];
vector[2] topleft_32[m[5]];
vector[2] bottomleft_32[m[5]];
//for three points case 4
vector[2] topright_33[m[6]];
vector[2] bottomright_33[m[6]];
vector[2] topleft_33[m[6]];
vector[2] bottomleft_33[m[6]];
// for 4 points
vector[2] topright[m[7]];
vector[2] bottomright[m[7]];
vector[2] topleft[m[7]];
vector[2] bottomleft[m[7]];
//for 2 points case 1
topright2[,1] = intervals_x_two[,2];
topright2[,2] = intervals_y_two[,2];
bottomright2[,1] = intervals_x_two[,2];
bottomright2[,2] = intervals_y_two[,1];
topleft2[,1] = intervals_x_two[,1];
topleft2[,2] = intervals_y_two[,2];
bottomleft2[,1] = intervals_x_two[,1];
bottomleft2[,2] = intervals_y_two[,1];
//for 2 points case 2
topright21[,1] = intervals_x_two1[,2];
topright21[,2] = intervals_y_two1[,2];
bottomright21[,1] = intervals_x_two1[,2];
bottomright21[,2] = intervals_y_two1[,1];
topleft21[,1] = intervals_x_two1[,1];
topleft21[,2] = intervals_y_two1[,2];
bottomleft21[,1] = intervals_x_two1[,1];
bottomleft21[,2] = intervals_y_two1[,1];
// for 3 points case1
topright_3[,1] = intervals_x_three[,2];
topright_3[,2] = intervals_y_three[,2];
bottomright_3[,1] = intervals_x_three[,2];
bottomright_3[,2] = intervals_y_three[,1];
topleft_3[,1] = intervals_x_three[,1];
topleft_3[,2] = intervals_y_three[,2];
bottomleft_3[,1] = intervals_x_three[,1];
bottomleft_3[,2] = intervals_y_three[,1];
//for 3 points case 2
topright_31[,1] = intervals_x_three1[,2];
topright_31[,2] = intervals_y_three1[,2];
bottomright_31[,1] = intervals_x_three1[,2];
bottomright_31[,2] = intervals_y_three1[,1];
topleft_31[,1] = intervals_x_three1[,1];
topleft_31[,2] = intervals_y_three1[,2];
bottomleft_31[,1] = intervals_x_three1[,1];
bottomleft_31[,2] = intervals_y_three1[,1];
//for 3 points case 3
topright_32[,1] = intervals_x_three2[,2];
topright_32[,2] = intervals_y_three2[,2];
bottomright_32[,1] = intervals_x_three2[,2];
bottomright_32[,2] = intervals_y_three2[,1];
topleft_32[,1] = intervals_x_three2[,1];
topleft_32[,2] = intervals_y_three2[,2];
bottomleft_32[,1] = intervals_x_three2[,1];
bottomleft_32[,2] = intervals_y_three2[,1];
//for 3 points case 4
topright_33[,1] = intervals_x_three3[,2];
topright_33[,2] = intervals_y_three3[,2];
bottomright_33[,1] = intervals_x_three3[,2];
bottomright_33[,2] = intervals_y_three3[,1];
topleft_33[,1] = intervals_x_three3[,1];
topleft_33[,2] = intervals_y_three3[,2];
bottomleft_33[,1] = intervals_x_three3[,1];
bottomleft_33[,2] = intervals_y_three3[,1];
//for 4 points
topright[,1] = intervals_x[,2];
topright[,2] = intervals_y[,2];
bottomright[,1] = intervals_x[,2];
bottomright[,2] = intervals_y[,1];
topleft[,1] = intervals_x[,1];
topleft[,2] = intervals_y[,2];
bottomleft[,1] = intervals_x[,1];
bottomleft[,2] = intervals_y[,1];
}
parameters{
vector[2] mu;
vector<lower=0>[2] sigma;
real<lower=-1, upper=1> rho;
}
transformed parameters{
cov_matrix[2] Sigma;
Sigma[1,1] = square(sigma[1]);
Sigma[1,2] = rho * sigma[1] * sigma[2];
Sigma[2,1] = Sigma[1,2];
Sigma[2,2] = square(sigma[2]);
}
model{
//for 2 points case 1
vector[m[1]] logg_2;//for cumulative CDF difference
//for 2 points case 2
vector[m[2]] logg_21;//for cumulative CDF difference
//for 2 pts
vector[2] topright_std_2[m[1]];
vector[2] bottomright_std_2[m[1]];
vector[2] topleft_std_2[m[1]];
vector[2] bottomleft_std_2[m[1]];
//for 2 pts case2
vector[2] topright_std_21[m[2]];
vector[2] bottomright_std_21[m[2]];
vector[2] topleft_std_21[m[2]];
vector[2] bottomleft_std_21[m[2]];
matrix[m[3], 4] con_mean3;
matrix[m[4], 4] con_mean31;
matrix[m[5], 4] con_mean32;
matrix[m[6], 4] con_mean33;
matrix[m[7], 4] con_mean;
vector[2] con_std;
//for 3 pts case1
vector[2] topright_std_3[m[3]];
vector[2] bottomright_std_3[m[3]];
vector[2] topleft_std_3[m[3]];
vector[2] bottomleft_std_3[m[3]];
//for 3 pts case2
vector[2] topright_std_31[m[4]];
vector[2] bottomright_std_31[m[4]];
vector[2] topleft_std_31[m[4]];
vector[2] bottomleft_std_31[m[4]];
//for 3 pts case3
vector[2] topright_std_32[m[5]];
vector[2] bottomright_std_32[m[5]];
vector[2] topleft_std_32[m[5]];
vector[2] bottomleft_std_32[m[5]];
//for 3 pts case4
vector[2] topright_std_33[m[6]];
vector[2] bottomright_std_33[m[6]];
vector[2] topleft_std_33[m[6]];
vector[2] bottomleft_std_33[m[6]];
//for 4 pts
vector[2] topright_std[m[7]];
vector[2] bottomright_std[m[7]];
vector[2] topleft_std[m[7]];
vector[2] bottomleft_std[m[7]];
//for 3 pts case 1
vector[m[3]] logg3;
//for 3 pts case 2
vector[m[4]] logg31;
//for 3 pts case 3
vector[m[5]] logg32;
//for 3 pts case 4
vector[m[6]] logg33;
//for 4 pts
vector[m[7]] logg;
//for 2 points case 1 and 2
for(i in 1:m[1]){
topright_std_2[i,1] = ( topright2[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topright_std_2[i,2] = ( topright2[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomright_std_2[i,1] = (bottomright2[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomright_std_2[i,2] = (bottomright2[i,2]- mu[2])/ sqrt(Sigma[2,2]);
topleft_std_2[i,1] = (topleft2[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topleft_std_2[i,2] = (topleft2[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomleft_std_2[i,1] = (bottomleft2[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomleft_std_2[i,2] = (bottomleft2[i,2]- mu[2])/ sqrt(Sigma[2,2]);
}
for(i in 1:m[2]){
topright_std_21[i,1] = ( topright21[i,2] - mu[1])/ sqrt(Sigma[1,1]);
topright_std_21[i,2] = ( topright21[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomright_std_21[i,1] = (bottomright21[i,2] - mu[1])/ sqrt(Sigma[1,1]);
bottomright_std_21[i,2] = (bottomright21[i,1]- mu[2])/ sqrt(Sigma[2,2]);
topleft_std_21[i,1] = (topleft21[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topleft_std_21[i,2] = (topleft21[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomleft_std_21[i,1] = (bottomleft21[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomleft_std_21[i,2] = (bottomleft21[i,1]- mu[2])/ sqrt(Sigma[2,2]);
}
//for 3 pts case 1
for(i in 1:m[3]){
topright_std_3[i,1] = (topright_3[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topright_std_3[i,2] = (topright_3[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomright_std_3[i,1] = (bottomright_3[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomright_std_3[i,2] = (bottomright_3[i,2]- mu[2])/ sqrt(Sigma[2,2]);
topleft_std_3[i,1] = (topleft_3[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topleft_std_3[i,2] = (topleft_3[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomleft_std_3[i,1] = (bottomleft_3[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomleft_std_3[i,2] = (bottomleft_3[i,2]- mu[2])/ sqrt(Sigma[2,2]);
}
//for 3 pts case 2
for(i in 1:m[4]){
topright_std_31[i,1] = (topright_31[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topright_std_31[i,2] = (topright_31[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomright_std_31[i,1] = (bottomright_31[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomright_std_31[i,2] = (bottomright_31[i,2]- mu[2])/ sqrt(Sigma[2,2]);
topleft_std_31[i,1] = (topleft_31[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topleft_std_31[i,2] = (topleft_31[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomleft_std_31[i,1] = (bottomleft_31[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomleft_std_31[i,2] = (bottomleft_31[i,2]- mu[2])/ sqrt(Sigma[2,2]);
}
//for 3 pts case 3
for(i in 1:m[5]){
topright_std_32[i,1] = (topright_32[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topright_std_32[i,2] = (topright_32[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomright_std_32[i,1] = (bottomright_32[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomright_std_32[i,2] = (bottomright_32[i,2]- mu[2])/ sqrt(Sigma[2,2]);
topleft_std_32[i,1] = (topleft_32[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topleft_std_32[i,2] = (topleft_32[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomleft_std_32[i,1] = (bottomleft_32[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomleft_std_32[i,2] = (bottomleft_32[i,2]- mu[2])/ sqrt(Sigma[2,2]);
}
//for 3 pts case 4
for(i in 1:m[6]){
topright_std_33[i,1] = (topright_33[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topright_std_33[i,2] = (topright_33[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomright_std_33[i,1] = (bottomright_33[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomright_std_33[i,2] = (bottomright_33[i,2]- mu[2])/ sqrt(Sigma[2,2]);
topleft_std_33[i,1] = (topleft_33[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topleft_std_33[i,2] = (topleft_33[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomleft_std_33[i,1] = (bottomleft_33[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomleft_std_33[i,2] = (bottomleft_33[i,2]- mu[2])/ sqrt(Sigma[2,2]);
}
//for 4 pts
for(i in 1:m[7]){
topright_std[i,1] = (topright[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topright_std[i,2] = (topright[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomright_std[i,1] = (bottomright[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomright_std[i,2] = (bottomright[i,2]- mu[2])/ sqrt(Sigma[2,2]);
topleft_std[i,1] = (topleft[i,1] - mu[1])/ sqrt(Sigma[1,1]);
topleft_std[i,2] = (topleft[i,2]- mu[2])/ sqrt(Sigma[2,2]);
bottomleft_std[i,1] = (bottomleft[i,1] - mu[1])/ sqrt(Sigma[1,1]);
bottomleft_std[i,2] = (bottomleft[i,2]- mu[2])/ sqrt(Sigma[2,2]);
}
//for 3 pts case 1
for(i in 1:m[3]){
con_mean3[i,1] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y_three[i,1] - mu[2]);//conditional mean for x given b1
con_mean3[i,2] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y_three[i,2] - mu[2]);//conditional mean for x given b2
con_mean3[i,3] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x_three[i,1] - mu[1]);//conditional mean for y given a1
con_mean3[i,4] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x_three[i,2] - mu[1]);//conditional mean for y given a2
}
//for 3 pts case 2
for(i in 1:m[4]){
con_mean31[i,1] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y_three1[i,1] - mu[2]);//conditional mean for x
con_mean31[i,2] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y_three1[i,2] - mu[2]);//conditional mean for x
con_mean31[i,3] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x_three1[i,1] - mu[1]);//conditional mean for y
con_mean31[i,4] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x_three1[i,2] - mu[1]);//conditional mean for y
}
//for 3 pts case 3
for(i in 1:m[5]){
con_mean32[i,1] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y_three2[i,1] - mu[2]);//conditional mean for x
con_mean32[i,2] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y_three2[i,2] - mu[2]);//conditional mean for x
con_mean32[i,3] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x_three2[i,1] - mu[1]);//conditional mean for y
con_mean32[i,4] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x_three2[i,2] - mu[1]);//conditional mean for y
}
//for 3 pts case 4
for(i in 1:m[6]){
con_mean33[i,1] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y_three3[i,1] - mu[2]);//conditional mean for x
con_mean33[i,2] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y_three3[i,2] - mu[2]);//conditional mean for x
con_mean33[i,3] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x_three3[i,1] - mu[1]);//conditional mean for y
con_mean33[i,4] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x_three3[i,2] - mu[1]);//conditional mean for y
}
for(i in 1:m[7]){
con_mean[i,1] = mu[1] + rho * sigma[1]/sigma[2] * (intervals_y[i,1] - mu[2]);//conditional mean for x given b1
con_mean[i,2] = mu[1] + rho * sigma[1]/sigma[2] *(intervals_y[i,2] - mu[2]);//conditional mean for x given b2
con_mean[i,3] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x[i,1] - mu[1]);//conditional mean for y given a1
con_mean[i,4] = mu[2] + rho * sigma[2]/sigma[1] *(intervals_x[i,2] - mu[1]);//conditional mean for y given a2
}
con_std[1] = sqrt(Sigma[1,1] * (1-square(rho)));
con_std[2] = sqrt(Sigma[2,2] * (1-square(rho)));
if(m[1]==0){
logg_2 = rep_vector(0, m[1]);
}
else{
for(i in 1:m[1]){
logg_2[i]= log(binormal_cdf(topright_std_2[i,1], topright_std_2[i,2], rho) - binormal_cdf(bottomright_std_2[i,1], bottomright_std_2[i,2], rho)
- binormal_cdf(topleft_std_2[i,1], topleft_std_2[i,2], rho) + binormal_cdf(bottomleft_std_2[i,1], bottomleft_std_2[i,2], rho)) ;
}
}
//(n-2) * log(binormal_cdf(topright_std_2[i,1], topright_std_2[i,2], rho) - binormal_cdf(bottomright_std_2[i,1], bottomright_std_2[i,2], rho)
//- binormal_cdf(topleft_std_2[i,1], topleft_std_2[i,2], rho) + binormal_cdf(bottomleft_std_2[i,1], bottomleft_std_2[i,2], rho))
//+ multi_normal_lpdf(bottomleft2[i,]| mu, Sigma) + multi_normal_lpdf(topright2[i,] | mu, Sigma);
//}
// to be added to the final log (n-2) * sum(logg_2) + multi_normal_lpdf(bottomleft2| mu, Sigma) + multi_normal_lpdf(topright2 | mu, Sigma);
if(m[2]==0){
logg_21 = rep_vector(0, m[2]);}
else{
for(i in 1:m[2]){
logg_21[i]= (n-2) * log(binormal_cdf(topright_std_21[i,1], topright_std_21[i,2], rho) - binormal_cdf(bottomright_std_21[i,1], bottomright_std_21[i,2], rho)
- binormal_cdf(topleft_std_21[i,1], topleft_std_21[i,2], rho) + binormal_cdf(bottomleft_std_21[i,1], bottomleft_std_21[i,2], rho))
+ multi_normal_lpdf(topleft21[i,]| mu, Sigma) + multi_normal_lpdf(bottomright21[i,] | mu, Sigma);
}
}
// to be added to the final log (n-2) * sum(logg_21) + multi_normal_lpdf(topleft21| mu, Sigma) + multi_normal_lpdf(bottomright21 | mu, Sigma);
//for 3 pts
if(m[3]==0){
logg3 = rep_vector(0, m[3]);}
else{
for(i in 1:m[3]){
logg3[i]= (n-3)* log(binormal_cdf(topright_std_3[i,1], topright_std_3[i,2], rho) - binormal_cdf(bottomright_std_3[i,1], bottomright_std_3[i,2], rho)
- binormal_cdf(topleft_std_3[i,1], topleft_std_3[i,2], rho) + binormal_cdf(bottomleft_std_3[i,1], bottomleft_std_3[i,2], rho))
+ multi_normal_lpdf(topright_3[i,]| mu, Sigma) + normal_lpdf(intervals_y_three[i,1]| mu[2], sigma[2])
+ normal_lpdf(intervals_x_three[i,1]| mu[1], sigma[1])
+(log(normal_cdf(intervals_x_three[i,2], con_mean3[i,1], con_std[1] ) - normal_cdf(intervals_x_three[i,1], con_mean3[i,1], con_std[1])) +
log(normal_cdf(intervals_y_three[i,2], con_mean3[i,3], con_std[2] ) - normal_cdf(intervals_y_three[i,1], con_mean3[i,3], con_std[2]))) ;
}
}
//for 3 pts version2
if(m[4]==0){
logg31 = rep_vector(0, m[4]);}
else{
for(i in 1:m[4]){
logg31[i]= (n-3)*log(binormal_cdf(topright_std_31[i,1], topright_std_31[i,2], rho) - binormal_cdf(bottomright_std_31[i,1], bottomright_std_31[i,2], rho)
- binormal_cdf(topleft_std_31[i,1], topleft_std_31[i,2], rho) + binormal_cdf(bottomleft_std_31[i,1], bottomleft_std_31[i,2], rho))
+multi_normal_lpdf(bottomleft_31[i,]| mu, Sigma)+ normal_lpdf(intervals_y_three1[i,2]| mu[2], sigma[2])
+ normal_lpdf(intervals_x_three1[i,2]| mu[1], sigma[1])+
(log(normal_cdf(intervals_x_three1[i,2], con_mean31[i,2], con_std[1] ) - normal_cdf(intervals_x_three1[i,1], con_mean31[i,2], con_std[1])) +log(normal_cdf(intervals_y_three1[i,2], con_mean31[i,4], con_std[2] ) - normal_cdf(intervals_y_three1[i,1], con_mean31[i,4], con_std[2])));
}
}
//for 3 pts case3
if(m[5]==0){
logg32 = rep_vector(0, m[5]);}
else{
for(i in 1:m[5]){
logg32[i]= (n-3)*log(binormal_cdf(topright_std_32[i,1], topright_std_32[i,2], rho) - binormal_cdf(bottomright_std_32[i,1], bottomright_std_32[i,2], rho)
- binormal_cdf(topleft_std_32[i,1], topleft_std_32[i,2], rho) + binormal_cdf(bottomleft_std_32[i,1], bottomleft_std_32[i,2], rho))
+ multi_normal_lpdf(bottomright_32| mu, Sigma) + normal_lpdf(intervals_y_three2[,2]| mu[2], sigma[2])
+ normal_lpdf(intervals_x_three2[,1]| mu[1], sigma[1]) +
(log(normal_cdf(intervals_x_three2[i,2], con_mean32[i,2], con_std[1] ) - normal_cdf(intervals_x_three2[i,1], con_mean32[i,2], con_std[1])) +log(normal_cdf(intervals_y_three2[i,2], con_mean32[i,3], con_std[2] ) - normal_cdf(intervals_y_three2[i,1], con_mean32[i,3], con_std[2])))
;
}
}
//for 3 pts case4
if(m[6]==0){
logg33 = rep_vector(0, m[6]);}
else{
for(i in 1:m[6]){
logg33[i]= (n-3)*log(binormal_cdf(topright_std_33[i,1], topright_std_33[i,2], rho) - binormal_cdf(bottomright_std_33[i,1], bottomright_std_33[i,2], rho)
- binormal_cdf(topleft_std_33[i,1], topleft_std_33[i,2], rho) + binormal_cdf(bottomleft_std_33[i,1], bottomleft_std_33[i,2], rho))
+ multi_normal_lpdf(topleft_33| mu, Sigma) + normal_lpdf(intervals_y_three3[,1]| mu[2], sigma[2])
+ normal_lpdf(intervals_x_three3[,2]| mu[1], sigma[1])+
(log(normal_cdf(intervals_x_three3[i,2], con_mean33[i,1], con_std[1] ) - normal_cdf(intervals_x_three3[i,1], con_mean33[i,1], con_std[1])) +log(normal_cdf(intervals_y_three3[i,2], con_mean33[i,4], con_std[2] ) - normal_cdf(intervals_y_three3[i,1], con_mean33[i,4], con_std[2])));
}
}
//for 4 pts
for(i in 1:m[7]){
logg[i]= (n-4) * log(binormal_cdf(topright_std[i,1], topright_std[i,2], rho) - binormal_cdf(bottomright_std[i,1], bottomright_std[i,2], rho)
- binormal_cdf(topleft_std[i,1], topleft_std[i,2], rho) + binormal_cdf(bottomleft_std[i,1], bottomleft_std[i,2], rho))
+
(log(normal_cdf(intervals_x[i,2], con_mean[i,1], con_std[1] ) - normal_cdf(intervals_x[i,1], con_mean[i,1], con_std[1])) +
log(normal_cdf(intervals_x[i,2], con_mean[i,2], con_std[1] ) - normal_cdf(intervals_x[i,1], con_mean[i,2], con_std[1])) +
log(normal_cdf(intervals_y[i,2], con_mean[i,3], con_std[2] ) - normal_cdf(intervals_y[i,1], con_mean[i,3], con_std[2])) +
log(normal_cdf(intervals_y[i,2], con_mean[i,4], con_std[2] ) - normal_cdf(intervals_y[i,1], con_mean[i,4], con_std[2])))+
(normal_lpdf(intervals_x[i,1]| mu[1], sigma[1])) +
(normal_lpdf(intervals_x[i,2]| mu[1], sigma[1])) +
(normal_lpdf(intervals_y[i,1]| mu[2], sigma[2])) +
(normal_lpdf(intervals_y[i,2]| mu[2], sigma[2]));
}
target +=
// for 2-point case1
(n-2) * sum(logg_2) + multi_normal_lpdf(bottomleft2| mu, Sigma) + multi_normal_lpdf(topright2 | mu, Sigma)
//for 2-point case2
+ sum(logg_21)
// for 3-point case1
+ sum(logg3)
// for 3-point case2
+ sum(logg31)
//for 3-point case3
+sum(logg32)
//for 3-point case3
+ sum(logg33)
//for 4 points
+ sum(logg) ;
}