Fitting data with Lmfit

I have two sets of data, one for the x axis and the other one for the y axis.

I'd like to make a fit of this data, with Lmfit, using the following definition of the Lèvy distribution.

equation

I have a constraint on c paramter, which means it must be greater than zero.

Below there's my code:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model
from lmfit import Parameters

def Levy(x, c, m):
    sq = np.sqrt(c/(2*np.pi))
    ex = np.exp(-c/(2*(x-m)))
    den = (x-m)**(3/2)
    return (sq*(ex/den))

x =np.array([0.03,0.08,0.13,0.18,0.23,0.28,0.33,0.38,0.43,0.48,0.53,0.58,0.63,0.68,0.73,0.78,0.83,0.88,0.93,0.98,1.03,1.08,1.13,1.18,1.23,1.28,1.33,1.38,1.43,1.48,1.53,1.58,1.63,1.68,1.73,1.78,1.83,1.88,1.93,1.98,2.03,2.08,2.13,2.18,2.23,2.28,2.33,2.38,2.43,2.48,2.53,2.58,2.63,2.68,2.73,2.78,2.83,2.88,2.93,2.98,3.03,3.08,3.13,3.28,3.88])

y = np.array([0.9931429,0.98486193,0.9783219,0.96919757,0.95680234,0.93760247,0.91618782,0.89456217,0.87109068,0.85210239,0.83358852,0.81407275,0.78769989,0.769608,0.73258027,0.71332794,0.68289386,0.65704847,0.62566436,0.60256181,0.57835128,0.5548792,0.53167057,0.50746062,0.48303911,0.45002014,0.4283945,0.40618835,0.38540666,0.36984661,0.34648062,0.3311846,0.30871501,0.29568682,0.28075974,0.2635118,0.25164404,0.23793043,0.22806706,0.21794024,0.2047541,0.1938894,0.18065023,0.16651464,0.15079665,0.14093328,0.1260062,0.1199406,0.11350606,0.1043281,0.09615262,0.08565687,0.07763934,0.07194268,0.0637672,0.05701618,0.04615091,0.03945234,0.03285927,0.03238484,0.0287456,0.0242624,0.01529543,0.00986279,0.00173977])

gmodel = Model(Levy)
params = Parameters()
params.add('c', value=0.2)
params.add('m', value=0)
result3 = gmodel.fit(y, x=x, params=params)

plt.plot(x,y,'bo')
plt.plot(x, result3.best_fit)
plt.show()

I have two problems:

How can I plot the fit for different vaues of c?
The peak of the Lèvy distribution does not reach the peak of the data set plot, even if I change the value of c.

Could someone help me, please? Thanks!

Solution

If I'm not misktaken, the µ parameter in the Lévy distribution is position parameter shifting the point where the probability density function starts to be non-null.

Given your data, if a Lévy distribution would describe it well then this parameter would be equal to 0 so we are left with one parameter, the scaling parameter c.

It seems that a Lévy distribution might not be the best probability distribution to describe your data. I rescaled your data so that the area under the curve would be equal to 1 (as for any probability density function) and plotted it together with the Lévy distribution with various values for c and indeed none seem to describe well your data.

f = plt.figure()
ax = plt.gca()
t = np.linspace(0, 4, 500)
for c in [0.2, 0.5, 1, 2, 4]:
    ax.plot(t, Levy(t, c, 0), label=r"$c=%s$"%c)

ax.set_xlabel("x")
ax.set_ylabel("y")

area_y = sum([(x[i+1]-x[i])*(y[i+1]+y[i])*0.5 for i in range(len(x)-1)])
y_ = y/area_y
ax.plot(x, y_, '.', label="data")

ax.legend()

If you want to fit such a function to data with constraints on the parameters you could alternatively use scipy optimize.curve_fit function.