用Lmfit拟合数据

Question

我有两组数据，一个是x轴，另一个是y轴。

我想用Lmfit对这些数据进行拟合，使用以下Lèvy发行版的定义。

​

​

我对c参数有一个限制，这意味着它必须大于零。

下面是我的代码：

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()

我有两个问题：

• 我怎样才能为不同的c值绘制合适的图呢？
• 即使更改c的值，Lèvy分布的峰值也不会达到数据集图的峰值。

有人能帮帮我吗？谢谢!

Answer 1

如果我没有错的话，Lévy分布中的偏差参数是位置参数，移动了概率密度函数开始为非零的点。

给定你的数据，如果Lévy分布能很好地描述它，那么这个参数将等于0，所以我们只剩下一个参数，缩放参数c。

似乎Lévy分布可能不是描述数据的最佳概率分布。我重新标度了你的数据，使曲线下的面积等于1(对于任何概率密度函数)，并将它与Lévy分布一起绘制，c的各种值，实际上，似乎没有一个能很好地描述你的数据。

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()

​

​

如果您想将这样的函数与具有参数约束的数据相匹配，您可以选择使用fit optimize.curve_fit函数。