Unexpected log_prob() behavior

[pgree]

I am getting unexpected behavior when I evaluate the log posterior using log_prob() in cmdstanpy. If I define the Stan program

parameters {
    real<lower=0> mu;
}
model {
    mu ~ normal(0, 1);
}

and then plot the posterior via

import numpy as np
import matplotlib.pyplot as plt
import cmdstanpy

# data of model
model = cmdstanpy.CmdStanModel(stan_file='model.stan')
n = 20
mus = np.linspace(0.0, 3.0, num=n)
post = np.zeros(n)
for i, mu in enumerate(mus):
    f1 = model.log_prob({'mu': mu})
    post[i] = np.exp(f1['lp_'])
plt.plot(mus, post, label='stan')
plt.legend()
plt.show()

then the resulting posterior is not normally distributed:

On the other hand, if I remove the constraint from the Stan program and instead use

parameters {
    real mu;
}
model {
    mu ~ normal(0, 1);
}

then I do get a normal distribution.

Any ideas why this might be happening? Thanks in advance for any help.

[jsocolar]

I don’t know for sure how cmdstanpy handles this stuff via log_prob, but it looks like you’re seeing the Jacobian adjustment here. There’s an awful lot of unconstrained parameter space packed into the area near mu = 0, so we need to adjust the target for this, which happens automatically when you declare the lower bound constraint.

[Bob_Carpenter]

@pgree

Any ideas why this might be happening?

Yes. The log density includes the Jacobian correction for the change of variables. See the Stan reference manual chapter on constrained variables and transforms. In this case of a positive constrained variable, the transform is mu = exp(mu_unconstrained), and we then multiply the density by the Jacobian of the transform, which is exp(mu_unconstrained). In other words, it’s like this:

parameters {
  real mu_unc;
}
transformed parameters {
  real<lower=0> mu = exp(mu_unc);
}
model {
  target += mu_unc;
  target += normal_lpdf(mu | 0, 1);
}

What you are getting out of log_prob is the log density on the unconstrained scale, which is the log density of mu_unc.

[pgree]

Thanks for your help @Bob_Carpenter and @jsocolar. Unfortunately, I’m still a bit confused. In this Stan program,

@Bob_Carpenter

parameters {
  real mu_unc;
}
transformed parameters {
  real<lower=0> mu = exp(mu_unc);
}
model {
  target += mu_unc;
  target += normal_lpdf(mu | 0, 1);
}

should mu_unc be distributed according to log_prob above, also the probability density function that I plotted here?

@pgree

then the resulting posterior is not normally distributed:

That appears to not be the case. Also, log_prob({'mu': mu}) can’t be evaluated at negative values of mu and the documentation, API Reference — CmdStanPy 1.2.0 documentation, says that the parameters inputted to log_prob should be given on the constrained scale.

[jsocolar]

When you evaluate model.log_prob({'mu': mu}), you are getting out the value of the log probability evaluated at the unconstrained parameter value that corresponds to the supplied value of mu (i.e. \mathcal{f}^{-1}(\mu), where \mathcal{f} is the constraining transform). The actual log probability density for small positive values of mu on the constrained scale is quite high. The reason this occurs is because of the distortion caused by the nonlinear constraining transform. Near zero, a small window on the constrained scale corresponds to a very large window on the unconstrained scale. Thus, narrow windows near zero on the constrained scale accumulate probability mass from a very broad window on the unconstrained scale. And thus, the probability density evaluated on the unconstrained scale at values that transform to near zero on the constrained scale is far less than the probability density evaluated on the constrained scale.

[pgree]

ok thanks @Bob_Carpenter @jsocolar

we want to sample X_1,\dots,X_n \sim X where X has p.d.f. f_X(x) for x > 0.

so we sample from Y_1, \dots, Y_n \sim Y = \log(X) and transform back via X_i = e^{Y_i} where Y has p.d.f.

f_Y(y) = f_X(g^{-1}(y)) | \frac{d}{dy} g^{-1}(y) | = f_X(e^y)e^y

for all y \in \mathbb{R} where g(y) = \log(y).

log_prob(x) is defined by \log(f_y(\log(x))). So if we want to evaluate the p.d.f. that we’re sampling from, we use the formula f_x(x) = f_Y(\log(x)) / x which can be evaluated via exp(log_prob(x)) / x