I’ll go through my through process so that there’s an example of how to track down fitting problems in models you haven’t used before.
Changing adapt_delta, and hence the accuracy of the symplectic integrator used by the Hamiltonian transition, didn’t affect the posterior fit which immediately suggested to me that there was a problem of numerical accuracy in your new density function instead of any geometry problems. Numerical inaccuracies lead to inaccurate gradients which break all of the wonderful properties of symplectic integrators, regardless of the step size; seeing the pathology unchanged as adapt_delta, and hence the step size, is typically a sign of numerical problems.
Indeed inspecting your normalhalflaplace_lpdf function one sees log1m(Phi(p)) which is an immediate red flag. For | p | \gtrapprox 5 the Phi function starts to saturate near 0 or 1 and loses precision which then causes the log1m evaluation to be significantly inaccurate. The resulting gradient is even worse. See https://github.com/stan-dev/math/issues/1284 for some recent discussion.
With that suspicion in mind I ran your code with data simulated from \alpha = 10 but printed out the (\alpha, \nu, \tau) configuration whenever | p | > 6. Every instance occurred for \alpha and \nu at their largest values, where you saw cutoffs in your posterior fits. In other words your posterior fit couldn’t go past those boundaries you saw because your log density function implementation wasn’t sufficiently accurate.
Then I went I a little bit too far and played around with asymptotic series to implement a more accurate log density function for p > 6. This ended up giving
real normalhalflaplace_lpdf(vector y, real alpha, real nu, real tau) {
int N = rows(y);
vector[N] lpd;
for(n in 1:N) {
real y_std = inv(tau) * (y[n] - nu);
real p = alpha * tau - y_std;
real log_Rp = 0;
if (p > 6) {
real inv_p = inv(p);
real inv_p_sq = square(inv_p);
real variable = 1;
real coefficient[7] = {-1, 2.5, -12.3333, 88.25, -816.2, 9200.83, -122028};
log_Rp = log(inv_p);
for (i in 1:7) {
variable *= inv_p_sq;
log_Rp += coefficient[i] * variable;
}
} else {
log_Rp = log1m(Phi(p)) - std_normal_lpdf(p);
}
lpd[n] = std_normal_lpdf(y_std) + log_Rp;
}
return sum(lpd) + N * log(alpha);
}
Running with this still gives a few divergences, but this time for small \alpha. The posterior, however, captures the simulated ground truth of \log(10) ~ 2.3.
The divergences for small \alpha could be due to numerical inaccuracies for p \lessapprox -5, but the pinching/funnel like shape suggests that it might also be inherent to the geometry of your problem and resolvable with a higher adapt_delta (the above plot was with default values). I also tried an asymptotic expansion for \log(1 - \Phi(x)) around x = \infty to account for possible numerical problems but couldn’t get anything to work.
Long story short: you were having fitting problems because of numerical instabilities in \log(1 - \Phi(x)) for x \gtrapprox 6. Adding a more robust implementation recovers a much better fit. That fit, however, demonstrates funnel-like behavior for small \alpha and small \nu, and poor identifiability between those two variables.
At the same time, however, the sampling of the data suggests another implementation. Simulating data with
real e1 = exponential_rng(1);
real z = normal_rng(0, 1);
y = nu + tau * z + e1 / alpha;
or mathematically,
suggests an expanded Stan program,
data {
int N;
vector[N] y;
}
parameters {
real<lower=0> alpha;
real epsilon[N];
real nu;
real<lower=0> tau;
}
model {
alpha ~ normal(0, 10);
nu ~ normal(0, 10);
tau ~ normal(0, 10);
epsilon ~ exponential(1);
y ~ normal(nu + epsilon / alpha, tau);
}
generated quantities {
real y_ppc[N] = normal_rng(nu + epsilon / alpha, tau);
}
In other words, don’t try to integrate the e1/\epsilon out of the model but fit everything jointly. That will result in a higher-dimensional (and more expensive) model but one without any serious numerical problems.