Modeling a Spectrum of Gaussians with SMC

[tvwenger]

@bob-carpenter

That’s one way to go, but it causes implicit label switching and thus can be hard to adapt and sample. An alternative that would be kosher according to the “generative” intentions of PyMC would be to do something like z ~ lognormal(mu, sigma) (or some other positive-constrained prior) and then set y = cumulative_sum(z). In Stan, we just transform with cumulative_sum(exp(x)) where x is unconstrained and apply the Jacobian correction and then put a prior directly on y.

Thank you for this suggestion, @bob-carpenter. I have experimented with such a model and have found good results!

coords = {"component": range(n_components)}
with pm.Model(coords=coords) as model:
    line_area_norm = pm.Gamma(
        "line_area_norm", alpha=2.0, beta=1.0, dims="component"
    )
    line_area = 100.0 * line_area_norm

    fwhm_norm = pm.Gamma(
        "fwhm_norm", alpha=2.0, beta=1.0, dims="component"
    )
    fwhm = 30.0 * fwhm_norm

    center_lower = -20.0
    center_offset_norm = pm.HalfNormal(
        "center_offset_norm",
        sigma=1.0,
        dims="component",
    )
    center_offset = 20.0 * center_offset_norm
    center = pm.Deterministic(
        "center",
        center_lower + pm.math.cumsum(center_offset, axis=0),
        dims="component"
    )
    
    noise = pm.HalfNormal("noise", sigma=1.0)

    amplitude = calc_amplitude(line_area, fwhm)
    predicted = gaussian(x[:, None], amplitude, fwhm, center).sum(axis=1)
    _ = pm.Normal(
        "predicted", mu=predicted, sigma=noise, observed=observed
    )

The model samples fine with both NUTS and SMC. Some NUTS chains collapse to a single model, but SMC seems to find both modes just fine.