python - 对数正态曲线拟合

Question

我在 x0 和 y0 数组中有对数正态分布数据:

x0.ravel() = array([19.8815 , 19.0141 , 18.1857 , 17.3943 , 16.6382 , 15.9158 ,
   15.2254 , 14.5657 , 13.9352 , 13.3325 , 12.7564 , 12.2056 ,
   11.679  , 11.1755 , 10.6941 , 10.2338 ,  9.79353,  9.37249,
    8.96979,  8.58462,  8.21619,  7.86376,  7.52662,  7.20409,
    6.89552,  6.6003 ,  6.31784,  6.04757,  5.78897,  5.54151,
    5.30472,  5.07812,  4.86127,  4.65375,  4.45514,  4.26506,
    4.08314,  3.90903,  3.74238,  3.58288,  3.4302 ,  3.28407,
    3.14419,  3.01029,  2.88212,  2.75943,  2.64198,  2.52955,
    2.42192,  2.31889,  2.22026,  2.12583,  2.03543,  1.94889,
    1.86604,  1.78671,  1.71077,  1.63807,  1.56845,  1.50181,
    1.43801,  1.37691,  1.31842,  1.26242,  1.2088 ,  1.15746,
    1.10832,  1.06126,  1.01619])

y0.ravel() =array([1.01567e+03, 8.18397e+02, 7.31992e+02, 1.11397e+03, 2.39987e+03,
   2.73762e+03, 4.65722e+03, 7.06308e+03, 9.67945e+03, 1.38983e+04,
   1.98178e+04, 1.97461e+04, 3.28070e+04, 4.48814e+04, 5.80853e+04,
   7.35511e+04, 8.94090e+04, 1.08274e+05, 1.28276e+05, 1.50281e+05,
   1.69258e+05, 1.91944e+05, 2.16416e+05, 2.37259e+05, 2.57426e+05,
   2.74818e+05, 2.90343e+05, 3.01369e+05, 3.09232e+05, 3.13713e+05,
   3.17225e+05, 3.19177e+05, 3.17471e+05, 3.14415e+05, 3.08396e+05,
   2.95692e+05, 2.76097e+05, 2.52075e+05, 2.29330e+05, 1.97843e+05,
   1.74262e+05, 1.46360e+05, 1.20599e+05, 9.82223e+04, 7.80995e+04,
   6.34618e+04, 4.77460e+04, 3.88737e+04, 3.23715e+04, 2.58129e+04,
   2.15724e+04, 1.58737e+04, 1.13006e+04, 7.64983e+03, 4.64590e+03,
   3.31463e+03, 2.40929e+03, 3.02183e+03, 1.47422e+03, 1.06046e+03,
   1.34875e+03, 8.26674e+02, 9.53167e+02, 6.47428e+02, 9.83651e+02,
   8.93673e+02, 1.23637e+03, 0.00000e+00, 8.36573e+01])

我想使用 curve_fit 来获取适合我的数据点的函数，以获得 mu(然后是中位数的 exp(mu))和该分布的 sigma。

import numpy as np
from scipy.optimize import *

def f(x, mu, sigma) :
   return 1/(np.sqrt(2*np.pi)*sigma*x)*np.exp(-((np.log(x)- 
   mu)**2)/(2*sigma**2))


params, extras = curve_fit(f, x0.ravel(), y0.ravel())

print "mu=%g, sigma=%g" % (params[0], params[1])

plt.plot(x0, y0, "o")
plt.plot(x0, f(x0 ,params[0], params[1])) 
plt.legend(["data", "fit"], loc="best")
plt.show()

结果如下:

mu=1.47897, sigma=0.0315236

Curve_fit

显然该函数无论如何都不适合数据。

当我将拟合函数乘以代码中的 1.3*10^(5) 时:

plt.plot(x0, 1.3*10**5*f(x0 ,params[0], params[1]))

这是结果:

manually changed fitting curve

计算出的 mu 值(相关正态分布的平均值)似乎是正确的，因为当我使用时:

 np.mean(np.log(x))

我得到1.4968838412183132，这与我从curve_fit获得的mu非常相似。

计算中位数

np.exp(np.mean(np.log(x))

给出 4.4677451525990675，这似乎没问题。

但是除非我看到拟合函数抛出了我的数据点，否则我并不真正相信这些数字。我的问题显然是，拟合函数没有(大) y0 值的信息。我怎样才能改变这一点？任何帮助表示感谢!

Answer 1

问题是，您的数据没有(!)显示对数正态 pdf，因为它们没有正确标准化。请注意，pdf 上的积分必须为 1。如果您对数据进行数值积分并以此进行标准化，例如

y1 = y0/np.trapz(x0, y0)

你的方法效果很好。

params, extras = curve_fit(f, x0, y1)

plt.plot(x0, y1, "o")
plt.plot(x0, f(x0 ,params[0], params[1])) 
plt.legend(["data", "fit"], loc="best")
plt.show()

和

print("mu=%g, sigma=%g" % (params[0], params[1]))

结果

mu=1.80045, sigma=0.372185