# Linear combination of two Gaussian functions: multimodal posterior

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*np.exp(-0.5*(x-x0)**2/w**2)
c2 = A*np.exp(-0.5*(x-x0)**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()
`````` 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)
``````

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.

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!

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