Linear combination of two Gaussian functions: multimodal posterior

Questions
1

Hi PyMC3 community,

Please forgive me it this is obvious, but I haven’t been able to find a clean workaround. I am new to Python and PyMC3.

I am running into a multimodal posterior distribution and want to restrict the parameter space so it becomes unimodal.

Below a MWE with an over-simplified model that describes a measured curve y=f(x) as a sum of two Gaussians over x, with fixed and identical width.

import pymc3 as pm
import numpy as np
import arviz as az
import matplotlib.pyplot as plt

# Define model: two Gaussians of identical fixed width over (0,1)
x = np.linspace(0,1,501)
w = 0.05
def y_model(x0,A):
    c1 = A[0]*np.exp(-0.5*(x-x0[0])**2/w**2)
    c2 = A[1]*np.exp(-0.5*(x-x0[1])**2/w**2)
    return c1+c2

# Simulate data
y = y_model([0.5,0.7],[0.7,0.3])
noise = np.random.normal(scale=0.02,size=np.size(x))
y_obs = y + noise

# Plot simulated data
plt.plot(x,y_obs)
plt.show()

image

Next the PyMC3 model and the MCMC simulation

# Probabilistic model with observed y
with pm.Model() as model:
    x0 = pm.Beta("x0",alpha=2,beta=2,shape=2)
    A = pm.Dirichlet("A",np.ones(2))
    sig = pm.HalfNormal("sig",sigma=0.1)
    y = pm.Normal("y",mu=y_model(x0,A),sigma=sig,observed=y_obs)

trace = pm.sample(model=model,cores=1)
pm.summary(trace)

image
image776×175 4.46 KB

The posterior mean of the positions x0 of the two Gaussian are equal and have a large spread; same for the amplitudes A. The R-hats are much larger than 1:

The trace plots reveal the problem: the posteriors for x0 and A are bimodal. That screws up the means, and the R-hats.

image
image994×499 101 KB

How to deal with this situation? How can I reformulate the PyMC3 model to render the posterior unimodal, for example by assigning index 0 always to the peak with smaller x0. The goal would be to get the posterior mean and spreads for each of the peaks.

Thanks!

2

I think this might answer your question: Fitting a spectra of gaussians