Hi all!
I‘m working on joint modeling of beta-binomial outcomes using Gauss Copula.
But stan warns that
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: Exception: inv_Phi: Probability variable is -2.14096e+06, but must be in the interval [0, 1] (in 'C:/Users/AAA/AppData/Local/Temp/RtmpKYJF6n/model-36d45da92db6.stan', line 19, column 6 to column 57) (in 'C:/Users/AAA/AppData/Local/Temp/RtmpKYJF6n/model-36d45da92db6.stan', line 129, column 1 to line 130, column 90)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine, Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
My model is:
y_1i~Beta-binomial(n_1i,alpha1,beta1);
y_2i~Beta-binomial(n_2i,alpha2,beta2);
(y1,y2)~Gauss_copula(L);
Here is my stan code：
The functions in the functions block except the beta_binomial_marginal could be found :here.
And the beta_binomial_marginal function is a modified version of binomial_marginal in the hyperlink.
functions{
//
array[] matrix beta_binomial_marginal(array[,] int num, array[,] int den, matrix alpha, matrix beta,
matrix u_raw) {
int N = rows(alpha);
int J = cols(alpha);
array[2] matrix[N, J] rtn;
for (j in 1 : J) {
for (n in 1 : N) {
real Ubound = beta_binomial_cdf(num[n, j] | den[n, j], alpha[n,j], beta[n,j]);
real Lbound = 0;
if (num[n, j] > 0) {
Lbound = beta_binomial_cdf(num[n, j]-1 | den[n, j], alpha[n,j], beta[n,j]);
}
real UmL = Ubound - Lbound; // Ubound>0,UmL>0
rtn[1][n, j] = inv_Phi(Lbound + UmL * u_raw[n, j]);
if(is_inf(rtn[1][n, j])) rtn[1][n, j] =negative_infinity();
if(is_nan(log(UmL))) rtn[2][n, j] = negative_infinity();
else rtn[2][n, j] = log(UmL);
}
}
return rtn;
}
real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) {
int N = rows(U);
int J = cols(U);
matrix[J, J] Gammainv = chol2inv(L);
return - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U))-N * sum(log(diagonal(L)));
}
real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) {
// Extract dimensions of final outcome matrix
int N = rows(marginals[1][1]);
int J = rows(L);
matrix[N, J] U;
// Iterate through marginal arrays, concatenating the outcome matrices by column
// and aggregating the log-likelihoods (from continuous marginals) and jacobian
// adjustments (from discrete marginals)
real adj = 0;
int pos = 1;
for (m in 1 : size(marginals)) {
int curr_cols = cols(marginals[m][1]);
U[ : , pos : (pos + curr_cols - 1)] = marginals[m][1];
adj += sum(marginals[m][2]);
pos += curr_cols;
}
// Return the sum of the log-probability for copula outcomes and jacobian adjustments
return multi_normal_cholesky_copula_lpdf(U | L) + adj;
}
}
data{
int<lower=0> N;
array[N,1] int<lower=0> y1;
array[N,1] int<lower=0> n1;
array[N,1] int<lower=0> y2;
array[N,1] int<lower=0> n2;
}
parameters{
real<lower=0,upper=1> a1; //a1=alpha1/(alpha1+beta1)
real<lower=1> b1; //b1=alpha1+beta1
real<lower=0,upper=1> a2;
real<lower=1> b2;
matrix<lower=0,upper=1>[N,1] u1;
matrix<lower=0,upper=1>[N,1] u2;
cholesky_factor_corr[2] L;
}
transformed parameters{
real<lower=0> alpha1;
real<lower=0> alpha2;
real<lower=0> beta1;
real<lower=0> beta2;
alpha1=a1*b1;
beta1=b1-alpha1;
alpha2=a2*b2;
beta2=b2-alpha2;
}
model{
//prior for a1,a2,b1,b2,rho
b1~pareto(1, 1.5);//
b2~pareto(1, 1.5);
L ~ lkj_corr_cholesky(1.0);
matrix[N,1] alpha_j1;
matrix[N,1] alpha_j2;
matrix[N,1] beta_j1;
matrix[N,1] beta_j2;
for(i in 1:N){
alpha_j1[i,1]=alpha1;
alpha_j2[i,1]=alpha2;
beta_j1[i,1]=beta1;
beta_j2[i,1]=beta2;
}
{beta_binomial_marginal(y1,n1,alpha_j1,beta_j1,u1),
beta_binomial_marginal(y2,n2,alpha_j2,beta_j2,u2)}~centered_gaussian_copula_cholesky(L);
//print("target=",target());
}
Here is my data and R code:
data<-read.table(text = "
37 50 28 44
49 65 33 55
59 79 46 64
25 33 24 21
59 79 45 64
13 18 12 14
25 34 17 31
27 35 23 34
107 141 99 108
76 101 64 78")
colnames(data)<-c("n1","n2","y1","y2")
stan_data<-list(N=dim(data)[1],
y1=as.matrix(data$y1),
n1=as.matrix(data$n1),
y2=as.matrix(data$y2),
n2=as.matrix(data$n2))
fit_betabinomial_gauss <- mod_betabinomial_gauss$sample(data = stan_data,#init = lapply(1:4,init_gen),
chains = 4,
parallel_chains = 4,
refresh = 0,
step_size = 0.01,adapt_delta=0.99,max_treedepth = 20)
Line 19 in the stan code is rtn[1][n, j] = inv_Phi(Lbound + UmL * u_raw[n, j]);
The Ubound and Lbound are calculated from the bulid-in function beta_binomial_cdf and u_raw(u1,u2) is bounded from 0 to 1 in the parameters block. It seems that Probability variable Lbound + UmL * u_raw[n, j] is bounded from 0 to 1 but in the warning it equals -2.14096e+06 .
Here is a screenshot in my testing:
If anyone has advice on how I can avoid getting this error, that’d be really appreciated.