scipy lognorm 不收敛于参数

Question

我已手动将对数正态分布拟合到我的数据中：

from scipy.stats import lognorm

sigma = 0.15
mu = 2
x_fit = np.linspace(x.min(), x.max(), 100)
y_fit = lognorm.pdf(x_fit, sigma, scale=np.exp(mu))
fig, ax = plt.subplots()
plt.plot(x_fit, y_fit, color='red')

问题是 scipy 不会收敛到 sigma = 0.15，但它只收敛到 mu ~ 2。另外，我也得到了一些 inf... 注意：即使我使用 [2, 0.15] 作为起始猜测，也会发生这种情况：

from scipy.optimize import curve_fit

y, x = np.histogram(z_data, bins=100)
x = (x[1:]+x[:-1])/2 # bin centers

def fit_lognormal(x, mu, sigma):
    return lognorm.pdf(x, sigma, scale=np.exp(mu))

p0 = [2, 0.15]
params = curve_fit(fit_lognormal, x, y, p0=p0)
params

返回：

(array([1.97713319e+00, 6.98238900e-05]),
 array([[inf, inf],
        [inf, inf]]))

Answer 1

假设您的数据集如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats, optimize

np.random.seed(123456)
sigma = 0.15
mu = 2

law = stats.lognorm(s=sigma, scale=np.exp(mu))
N = 30_000
x = law.rvs(size=N)

density, bins = np.histogram(x, density=1., bins=30)  # Mind the density=1.
centers = (bins[1:] + bins[:-1]) / 2.

您可以通过三种简单的方法来拟合参数：

原生

stats

:

stats.lognorm.fit(x)
(0.151297129052073, 0.053788977866467386, 7.335249847775119)

哪里

loc ~ 0

和

scale ~ np.exp(2)

。

使用

curve_fit

来适应 PDF：

def model(x, mu, sigma):
    return stats.lognorm.pdf(x, s=sigma, scale=np.exp(mu))

popt, pcov = optimize.curve_fit(model, centers, density)
# (array([2.00059957, 0.15158351]),
#  array([[8.87310596e-07, 4.50003479e-08],
#        [4.50003479e-08, 5.93772296e-07]])) 

使用 MLE 和

minimize

:

def likelihood(p):
    return - np.sum(stats.lognorm.logpdf(x, s=p[1], scale=np.exp(p[0])))

sol = optimize.minimize(likelihood, x0=[1, 1])
#     fun: 45693.61136284961
#  hess_inv: array([[ 3.38213795e-08, -1.68343490e-08],
#        [-1.68343490e-08,  5.54538464e-07]])
#       jac: array([-0.00048828,  0.00097656])
#   message: 'Desired error not necessarily achieved due to precision loss.'
#      nfev: 66
#       nit: 17
#      njev: 22
#    status: 2
#   success: False
#         x: array([2.00008072, 0.15018327])