2
$\begingroup$

In this PyMC3 tutorial on Bayesian Mixed Effects Models, there is some Re-parameterization "to avoid chain divergences."

with pm.Model(coords=coords) as hierarchical:
    # Hyperpriors
    intercept_mu = pm.Normal("intercept_mu", 0, sigma=1)
    intercept_sigma = pm.HalfCauchy("intercept_sigma", beta=2)
    slope_mu = pm.Normal("slope_mu", 0, sigma=1)
    slope_sigma = pm.HalfCauchy("slope_sigma", beta=2)

    # Define priors
    sigma = pm.HalfCauchy("sigma", beta=2, dims="group")

    # Reparameterise to avoid chain divergences
    # β0 = pm.Normal("β0", intercept_mu, sigma=intercept_sigma, dims="group")
    β0_offset = pm.Normal("β0_offset", 0, sigma=1, dims="group")
    β0 = pm.Deterministic("β0", intercept_mu + β0_offset * intercept_sigma, dims="group")

    # Reparameterise to avoid chain divergences
    # β1 = pm.Normal("β1", slope_mu, sigma=slope_sigma, dims="group")
    β1_offset = pm.Normal("β1_offset", 0, sigma=1, dims="group")
    β1 = pm.Deterministic("β1", slope_mu + β1_offset * slope_sigma, dims="group")

    # Data
    x = pm.Data("x", data.x, dims="observation")
    g = pm.Data("g", data.group_idx, dims="observation")
    # Linear model
    μ = pm.Deterministic("μ", β0[g] + β1[g] * x, dims="observation")
    # Define likelihood
    pm.Normal("y", mu=μ, sigma=sigma[g], observed=data.y, dims="observation")

Consider that β1 = pm.Normal("β1", slope_mu, sigma=slope_sigma, dims="group") was swapped out for the below:

β1_offset = pm.Normal("β1_offset", 0, sigma=1, dims="group")
    β1 = pm.Deterministic("β1", slope_mu + β1_offset * slope_sigma, dims="group")

Why is this necessary and what does it accomplish?

Edit: I see the comment link to Stan forum, is there a slightly more concise/higher level explanation available?

I understand that HMC doesn't do well with high curvature areas; but why is multiplying the slope by an additional Normally distributed random variable helpful?

$\endgroup$
2

1 Answer 1

1
+50
$\begingroup$

The explanation below is based on Statistical Rethinking, 2nd edition, pages 420-423.


As you well say, samplers have issues exploring high curvature areas. A re-parametrization is a way to change the area that is being explored.

Consider a funnel made from the following parametrization:

$$ \begin{aligned} v &\sim \operatorname{Normal}(0, 3) \\ x &\sim \operatorname{Normal}(0, \exp(v)) \end{aligned} $$

Where does the funnel come from? Half the values of $v$ will be negative, and for negative numbers $\exp(v) < 1$. So we know that $x$ will be centered around 0, with a quite small variance. For half the other values of $v$, $\exp(v)$ will not only be larger than 1, it will grow—literally—exponentially. Hence, a funnel arises: increasingly more concentrated values of $x$ for more negative values of $y$, increasingly more spread values of $x$ for more positive values of $y$

Consider a reparametrization now, as follows:

$$ \begin{aligned} v &\sim \operatorname{Normal}(0, 3) \\ z &\sim \operatorname{Normal}(0, 1) \\ x &= z \exp(v) \end{aligned} $$

Where did the funnel go? We are no longer sampling $x$, but a standard normal variable $z$. Gladly, multiplying a normally distributed variable by a constant, changes its scale (or variance). Hence we can get our desired $x$ by multiplying $z$ by the exponentiated $v$. This step is a deterministic one, there is no probability involved, so there is no need for the sampler to explore anything here. We see that the sampler now just has to explore normal distributions, which have straightforward curvatures. The joint distribution of both $v$ and $z$ will now just be a multivariate normal, which has nothing pointy or weird going on.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.