I have written quite a simple model in which the likelihood is a function of three parameters \mu, \sigma, \xi that are derived by a 1-1 continuous transformation of three independent parameters q_{1}, q_{2}, q_{3} each of which is given a gamma prior. I have written a custom probability function.

When I have tested this Stan program with 400 iterations and one chain, it occasionally converges very quickly to what I know by other means is definitely the right result, but mostly the sampler quickly locks on to arbitrary values and stays there, reporting that maximum tree depth has been reached.

I have just tried again with four chains. One of them mixed well and gave results that I know are reasonable. The other three did not. I looked at the summary sampler parameters in shinystan. The good chain shows:

accept_stat = 0.9989

stepsize = 0.0159

treedepth = 6.725

n_leapfrog = 206.5

energy = 5.06

By contrast, for example, the figures for one of the other chains were:

accept_stat = 0.11

stepsize = 0.0000

treedepth = 5.1

n_leapfrog = 53209

energy = -366

All three bad chains reported numerous cases of maximum treedepth being reached (I had set it to 20).

The Stan program is:

```
functions {
vector lambda(vector y, real mu , real sigma,real xi){
vector[rows(y)] result;
real comp;
if(fabs(xi) > 1e-15){
real inv_xi = inv(xi);
for(i in 1: rows(y)){
comp = 1 + xi * (y[i] - mu) / sigma;
result[i] = comp > 0 ? -(1 + inv_xi) * log(comp) - log(sigma): 0.0;
}
}
else{
for(i in 1: rows(y)){
comp = (y[i] - mu) / sigma;
result[i] = comp > 0 ? comp - log(sigma): 0.0;
}
}
return result ;
}
real Blambda(
real mu,
real sigma,
real xi,
real nyears,
real u
){
real inv_xi = inv(xi);
real term1;
if (fabs(xi) > 1e-15){
term1 = 1 + xi*(u - mu) / sigma;
}
else{
term1 = exp((u - mu) / sigma);
}
term1 = term1 > 0 ? - ((term1) ^ (- inv_xi) )* nyears: 0;
return term1;
}
real poispt_lpdf(
vector y,
real mu,
real sigma,
real xi,
real nyears,
real u
){
real L;
real l;
L = Blambda(mu, sigma, xi, nyears, u);
l = sum(lambda(y, mu, sigma, xi));
return L + l;
}
}
data {
int<lower=0> N;
vector[N] y;
real nyears;
real u;
}
parameters {
real<lower = 0> qtilde1;
real<lower = 0> qtilde2;
real<lower = 0> qtilde3;
}
transformed parameters{
real mu;
real xi;
real<lower = 0> sigma;
real<lower = 0> nu;
nu = qtilde3 / qtilde2;
xi = log10(nu);
sigma = qtilde2 * xi / (nu * (nu - 1));
mu = qtilde1 - qtilde2 / nu;
}
model {
// priors
qtilde1 ~ gamma(38.9, 0.67);
qtilde2 ~ gamma(7.1, 0.16);
qtilde3 ~ gamma(47.0, 0.39);
// includes Jacobian adjustment
target += poispt_lpdf(y | mu, sigma, xi, nyears, u) -
log(xi) + log(nu) + log(nu - 1) + log(qtilde3);
}
```