Sampler appears to avoid zero in truncated hierarchical priors

Summary:

Samples generated from a simple hierarchical normal prior (no model, just prior) don't agree with the theoretical distribution or those generated numerically near the truncation point (zero)

Description:

This is from a Discourse thread: http://discourse.mc-stan.org/t/sampler-avoiding-zero-in-hierarchical-prior/

I’m comparing the results from using hmc-nuts to sample directly from a (hierarchical) prior to using normal_rng to generate samples. What I’ve found is that for the scale parameter of the normal distribution, the sampled prior doesn’t agree with the generated one - and for some reason the sampled prior avoids zero.

Here's the code for a minimal example:

parameters {
  real<lower=0> sigma_bar;
  real<lower=0> sigma_sig;
  real<lower=0> sigma;  
}

model {
  sigma_bar ~ normal(0.5,0.5) T[0,];
  sigma_sig ~ normal(0.5,0.5) T[0,];
  sigma ~ normal(sigma_bar,sigma_sig) T[0,]; 
}

generated quantities {
 real<lower=0> s_bar;
 real<lower=0> s_sig;
 real<lower=0> s;
s_bar = -1;
s_sig = -1;
s = -1;
while(s_bar <0)
 s_bar = normal_rng(0.5,0.5);
while(s_sig < 0)
 s_sig = normal_rng(0.5,0.5);
while(s < 0)
 s = normal_rng(s_bar,s_sig);
}

Plotting the output, I find that the distribution of s_sig (the generated quantity) matches the theoretical PDF for the truncated normal perfectly, but sigma_sig (the sampled quantity) has a discrepancy near zero.

normal

Reproducible Steps:

I'm running MatlabStan - this is my code to generate that figure, more or less - stanfile.stan is the code listed above. Note to be super confident in these histograms, I generate a bunch of samples (2e7) - this takes about ten minutes to run on my machine, creates csv files close to a gig.

dat = struct;
fit = stan('file','stanfile.stan','data',dat,'iter',2e7,'chains',4);
fit.verbose=true
  fit.block()
z = fit.extract('pars',{'sigma_bar','sigma_sig','sigma','s','s_bar','s_sig'});
[n1,x1] = hist(z.s_sig,100);
[n2,x2] = hist(z.sigma_sig,100);
plot(x1,n1/trapz(x1,n1),'LineWidth',2);
hold on
plot(x2,n2/trapz(x2,n2),'LineWidth',3,'color',[0.7 0.7 0.7]);
phi = @(z) (1/sqrt(2*pi))*exp(-0.5*z.^2);
PHI = @(x) (1/2)*(1+erf(x/sqrt(2)));
mu=0.5; sigma=0.5; % params for the truncated normal
pdf_trunc = phi((x2-mu)/sigma)/(sigma*(1-PHI((-mu/sigma))));
plot(x2,pdf_trunc,'r--','LineWidth',2);
set(gca,'FontSize',18);
xlabel('sigma_sig','interpreter','none');
ylabel('Probability density')
legend({'Generated (s_sig)','Sampled (sigma_sig)','True PDF'},'interpreter','none');

Current Version:

I'm doing this with v2.14.0 (this is post sampler-bug-fix) and I don't think the sampler has changed since then. It shouldn't matter, hopefully, but I'm also using MatlabStan

Start with the minimal hierarchical Stan model in the issue and reproduce it with MatlabStan v2.14.0, comparing sigma_sig with the generated s_sig and the stated truncated-normal PDF. Read the linked Discourse thread for the existing investigation; done means explaining the discrepancy near zero and establishing whether the sampler or comparison is responsible.