I have written my own _lpdf function, defined it in the functions block, and it is at least syntactically correct. The function is rather complex, as I had to deal with ragged data, though I believe that to be irrelevant for my question.

Initializing with random values or 0, I get:

`Log probability evaluates to log(0)`

So I added

`print(lprob)`

to my code, where lprob is a real scalar which is returned by the function.

I copied the console output into a text editor, looking for ‘nan’, ‘na’, and also inspected the output just by scrolling through the numbers. I cannot find a single nan. All numbers look fine to me.

What else can I do to debug?

Full code, sorry for all the idiosyncrasies. The unit/household/person/i-specific log-likelihood contains PDF, CDF components + Jacobian from COV arithmetic.

Standard multivariate Normal heterogeneity on top of it.

```
functions{
real vdm_lpdf(vector xx, vector xxp, vector pp, vector aai,
real mxp, vector minz,
vector betai, real gamma, real E, real sigma,
int ntski, int nalti, int nattr, int dimxi){
//defs
vector[dimxi] ab; // where a*beta is saved
vector[dimxi] g; // store g vector
vector[dimxi] j1; // Jacobian part 1
real j2; // Jacobian part 2
vector[dimxi] j21;
vector[nalti] j22;
real lprob; //variable to hold log prob
int pos; // counters
int pos2;
int pos3;
//compute a*beta
pos=1;
pos2=1;
for(jj in 1:dimxi){
ab[pos2]=sum(betai .* segment(aai,pos,nattr));
pos=pos+nattr;
pos2=pos2+1;
}
//g vector
g = -ab + log(gamma*xx+ 1) + log(pp) - log(minz+E-mxp) + log(pp);
//Jacobian
j1 = (log(gamma)-log( (gamma * xx) + 1)) .* (xxp);
j21 = ( ((gamma*xx+1) .* pp) ./ (gamma*(minz+E-mxp))) .* (xxp);
pos3=1;
for(jj in 1:ntski){
j22[jj] = sum(segment(j21,pos3,nalti));
pos3=pos3+nalti;
}
j2=sum(log(j22+1));
//ll
lprob = -sum(exp(-g/sigma)) + sum((g/sigma-log(sigma)) .* xxp) + sum(j1) + j2;
//print(lprob);
return lprob;
}
}
data{
int<lower=1> N; // no of households
int<lower=1> K; // no of hh level parameters
int<lower=1> nattr; // no of attributes
int<lower=1> dimxst; // total number of elements for x, p
int<lower=1> dimast; // total number of elements for a
int<lower=1> dimxs[N]; // hh's number of elements for x
int<lower=1> dimas[N]; // hh's number of elements for a
int<lower=1> ntsks[N]; // number of choice tasks for household
int<lower=1> nalts[N]; // number of alternatives available to household
vector<lower=0>[N] mxps; // max spent for hh
vector<lower=0>[dimxst] XX; // all choices X in one vector
vector<lower=0>[dimxst] XXP; // all choices X in one vector
vector<lower=0>[dimxst] PP; // all prices in one vector
vector<lower=0>[dimxst] MINZ; // minimum outside good z for each task/alternative for each hh
vector[dimast] AAI; // all attributes in one vector
}
parameters{
matrix[nattr,N] beta;
vector[N] lgammas;
vector[N] lsigmas;
vector<lower=1>[N] Eus;
corr_matrix[K] Omega;
vector<lower=0>[K] tau;
vector[K] MU;
}
transformed parameters{
vector[N] SIGMAS;
vector[N] GAMMAS;
cov_matrix[K] SIGMA;
vector[N] Es;
matrix[K,N] theta;
vector[N] LE;
LE = Eus .* log(mxps);
Es = exp(LE)+0.00001;
theta = append_row(append_row(append_row(beta,lgammas'),LE'),lsigmas');
SIGMA = quad_form_diag(Omega, tau);
SIGMAS = exp(lsigmas);
GAMMAS = exp(lgammas);
}
model{
//counters
int posm;
int posm2;
posm=1;
posm2=1;
//priors
tau ~ cauchy(0, 2.5);
Omega ~ lkj_corr(2);
MU ~ normal(0,10);
for (j in 1:N){
theta[,j] ~ multi_normal(MU, SIGMA);
}
//hh level lls
for(i in 1:N){
segment(XX,posm,dimxs[i]) ~ vdm(segment(XXP,posm,dimxs[i]),
segment(PP,posm,dimxs[i]),
segment(AAI,posm2,dimas[i]),
mxps[i],
segment(MINZ,posm,dimxs[i]),
beta[,i],
GAMMAS[i],
Es[i],
SIGMAS[i],
ntsks[i],
nalts[i],
nattr,
dimxs[i]);
posm=posm+dimxs[i];
posm2=posm2+dimas[i];
print(target());
}
}
```