Simpel model to estimate loaction of Lorentzian peak not working

Questions v3
1

Hi,

I am building a simple model to infer the location of a (narrow) Lorentzian shaped peak and run into problems. To illustrate the issues I first generate some data using a Lorentzian peak (line) model:

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0.0, 2000.0, 4096)
sigma_noise = 0.002 # sig for gaussian noise
width = 1.0         # with of peak
x_target = 400.0    # true location of the peak
a0 = 3

# Lorentzian line at x0, int_{-oo}^{oo} line(x,x0)dx = a0
def line(x, x0):
    return 2*a0/(1+(2*(x-x0)/width)**2)/(np.pi*width)

# Generate some data
data = line(x, x_target)
# add gaussian noise
noise = np.random.normal(0, sigma_noise, size=x.size)
data += noise

# Plot data
plt.figure(1)
plt.plot(x, data, '-k', label='data')
plt.plot(x, line(x,440), '-b', label='model peak at x0=440')
plt.xlabel('x')
plt.legend()

data

The typical uncertainty in the predicted peak location x_pred is around 100 times the width of the peak. (Note that the x range extents far beyond that. This is because later in a more complex model, there could be other peaks at any location.)

When I try to estimate the location of the peak using the simple PyMC3 model below:

import arviz as az
import pymc3 as pm

print(f"Running on PyMC3 v{pm.__version__}")
print(f"Running on ArviZ v{az.__version__}")

#x_pred_mu = x_target
x_pred_mu = 440
x_pred_err = 80

with pm.Model():

    x_pred = pm.Normal('x_pred', mu=x_pred_mu, sigma=x_pred_err)    
    
    # residual
    res = line(x, x_pred) - data
    pm.Normal("likelihood", mu=res, sigma=sigma_noise, observed=0)

    trace = pm.sample(draws=2000, tune=1000)

    az.style.use("arviz-darkgrid")    
    az.plot_trace(trace, var_names=["x_pred"], lines=[("x_pred", {}, x_target)])

Sampling fails producing the following results:

res1
res11187×509 82.2 KB

Using a truncated prior distribution, I have tried a TruncatedNormal and a Uniform distribution, did not produce better results. Any help or input would be greatly appreciated.