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