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
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.