R encounter fatal error and session terminated while using Rstan

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

How are you calling stan?

I first compiled the model
sym_bivariate_mle <-stan_model(file=“bivariate_interval_minmax_combo1.stan”)

and then use the optimisation function
sym_mle_correct_data_combo[[i]] <-optimizing(sym_bivariate_mle, data=data_set1,init=initf1)$par

calling stan, I did library(rstan)

I have tried the same 45 combo on three other likelihood formulations and the optimisation runs okay for all three likelihood construction, I suspect either the above likelihood function causes the termination of R ( but why it works for some cases? ) or there is something I overlook?

It is not supposed to crash R ever, although that happens sometimes. What is the 20th scenario?

20th scenario just an example, it definitely runs on the last scenario where N (the number of classical datapoints are 100), Rho = 0.9 and M( number of “symbols” , bivariate interval or rectangular representing min max of these classical data are 50).

I tried it again, it runs on the 45th, but it didn’t run on the 1st one which N =60, Rho= 0 and M =20

I am experiencing the same problem - the Stan model seems to be fine but the R session crashes early on during sampling (chain 1) with the warning “R session encountered a fatal error. R session was terminated” . Any ideas why this might be happening? I run the model in Rstan on Windows, stan model is attached. Any help is much appreciated !! Laura

I call the model in R:
fit <- stan(file = mname, data = inp, iter = niter, warmup = nwarmup, chains = nchains, verbose = TRUE, refresh = nrefresh)

revL.stan (2.3 KB)

Assuming you are using R through Rstudio does the behavior differ you set options(mc.cores = parallel::detectCores()) or not?