This problem has me completely stumped. I wrote a likelihood based on a Gaussian copula with beta marginals. For now, I’m working on a 2-dimensional copula, so I only need a matrix of means (2 for each observation) mu, a vector of precisions phi, and a correlation parameter rho. This is the function:
real gcbm_lpdf (matrix y, matrix mu, vector phi, real rho) {
// Variable declarations
int N = rows(y);
int D = cols(y);
real a;
real b;
vector[D] z;
vector[D] ldy;
// Log-likelihood
real r2 = square(rho);
real ldr = log1m(r2);
real rr2 = r2/(1 - r2);
vector[N] ll;
for (i in 1:N) {
for (j in 1:D) {
a = mu[i, j]*phi[j];
b = (1 - mu[i, j])*phi[j];
z[j] = inv_Phi(inc_beta(a, b, y[i, j]));
ldy[j] = (a - 1)*log(y[i, j]) + (b - 1)*log1m(y[i, j]) - lbeta(a, b);
}
ll[i] = -0.5*(ldr + rr2*(square(z[1]) - 2*z[1]*z[2]/rho + square(z[2]))) + sum(ldy);
}
return sum(ll);
}
The means I obtain for a particular model and since they must lie between 0 and 1, I use the following (an almost ideal demand system):
matrix aid (vector e, matrix p, real alpha0, row_vector alpha, matrix gamma, row_vector[] beta, int trm) {
// Variable declarations
int N = rows(p);
int d = cols(p);
int D = d - 1;
real elap;
matrix[trm, d] sVec;
matrix[N, d] mu;
// Mean computation
matrix[N, d] pg = p*gamma;
for (i in 1:N) {
elap = e[i] - alpha0 - dot_product(p[i], alpha) - 0.5*quad_form(gamma, p[i]');
for (r in 1:trm) {
sVec[r] = beta[r]*pow(elap, r);
}
mu[i] = alpha + pg[i] + rep_row_vector(1, trm)*sVec;
}
return mu[, :D];
}
My big issue is that I get the following error every time I try to sample from this model:
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: Exception: Error in function boost::math::ibeta<long double>(long double, long double, long double): The argument a to the incomplete beta function must be >= zero (got a=-185.490376487676513761). (in 'Code/Stan/Functions.stan' at line 71; included from 'model29609967fdf_Gaussian_Beta_AID' at line 1)
(in 'model29609967fdf_Gaussian_Beta_AID' at line 47)
However, I have checked using R, expose_stan_functions(), and the generated quantities block that I get a valid likelihood using exactly this function and the initial parameter values I supply. Does anyone have any insight how might I debug this issue? I’ve tried everything and it is the first time I’m running into this problem after successfully using Stan several times.
Edit: I should have added that since the parameters are either centered or unit-sum vectors, I supply initial values for the raw parameters and then obtain the full parameters in the transformed parameters block. Here is the full model just in case (right now the priors are extremely tight because I thought that might help but that will change in the actual estimation):
data {
int N;
int D;
int d;
int dg;
int trm;
matrix[N, D] y;
vector[N] e;
matrix[N, d] p;
}
parameters {
real<lower=0> a0;
row_vector<lower=0>[D] a;
row_vector<upper=0>[D] b[trm];
vector[dg] g;
vector<lower=0>[D] phi;
real<lower=-1, upper=1> rho;
}
transformed parameters{
row_vector[d] alpha;
matrix[d, d] gamma;
row_vector[d] beta[trm];
alpha = append_col(a, 1 - sum(a));
gamma = g_to_gamma(g); // Transform parameter vector to symmetric matrix
for (r in 1:trm) {
beta[r] = append_col(b[r], 0-sum(b[r]));
}
}
model {
matrix[N, D] mu;
// Priors
a0 ~ normal(1, 0.1);
a[1] ~ normal(0.25, 0.01);
a[2] ~ normal(0.25, 0.01);
g ~ normal(0, 0.01);
for (r in 1:trm) {
to_vector(b[r]) ~ normal(0, 0.01);
}
phi ~ gamma(10, 1);
rho ~ uniform(-1, 1);
// Likelihood
mu = aid(e, p, a0, alpha, gamma, beta, trm);
y ~ gcbm(mu, phi, rho);
}