Problem with rejection of initial values in custom likelihood

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

Are those warnings from warmup period or from sampling period?

Do you give inits for the sampling call?

They are from before the warmup period starts, as I cannot get the sampling started due to an error in the initial values. I am passing them through inits via the function method, as I perturbate some initial values that I have found give a finite likelihood. However, even if I pass the regular values, it still does prints that same error. I believe they must be related to the values of what I call b in the model. In the expose_stan_functions() it works only when I pass them as lists of row vectors, so this is how I pass them to inits.