使用 CURVE_FIT 在 Python 中拟合对数正态分布

Question

我有一个 x 的假设 y 函数，并试图找到/拟合一条对数正态分布曲线，该曲线可以最好地塑造数据。我正在使用 curve_fit 函数并且能够拟合正态分布，但曲线看起来并不优化。

下面是给出 y = f(x) 的 y 和 x 数据点。

y_axis = [0.00032425299473065838, 0.00063714106162861229, 0.00027009331177605913, 0.00096672396877715144, 0.002388766809835889, 0.0042233337680543182, 0.0053072824980722137, 0.0061291327849408699, 0.0064555344006149871, 0.0065601228278316746, 0.0052574034010282218, 0.0057924488798939255, 0.0048154093097913355, 0.0048619350036057446, 0.0048154093097913355, 0.0045114840997070331, 0.0034906838696562147, 0.0040069911024866456, 0.0027766995669134334, 0.0016595801819374015, 0.0012182145074882836, 0.00098231827111984341, 0.00098231827111984363, 0.0012863691645616997, 0.0012395921040321833, 0.00093554121059032721, 0.0012629806342969417, 0.0010057068013846018, 0.0006081017868837127, 0.00032743942370661445, 4.6777060529516312e-05, 7.0165590794274467e-05, 7.0165590794274467e-05, 4.6777060529516745e-05]

y 轴是事件在 x 轴时间段中发生的概率：

x_axis = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0]

我能够使用 excel 和对数正态方法更好地拟合我的数据。当我尝试在 python 中使用对数正态时，拟合不起作用，我做错了什么。

以下是我用于拟合正态分布的代码，这似乎是我可以在 python 中拟合的唯一代码（难以置信）：

#fitting distributino on top of savitzky-golay
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import scipy
import scipy.stats
import numpy as np
from scipy.stats import gamma, lognorm, halflogistic, foldcauchy
from scipy.optimize import curve_fit

matplotlib.rcParams['figure.figsize'] = (16.0, 12.0)
matplotlib.style.use('ggplot')
# results from savgol
x_axis = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,     13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0]
y_axis = [0.00032425299473065838, 0.00063714106162861229, 0.00027009331177605913, 0.00096672396877715144, 0.002388766809835889, 0.0042233337680543182, 0.0053072824980722137, 0.0061291327849408699, 0.0064555344006149871, 0.0065601228278316746, 0.0052574034010282218, 0.0057924488798939255, 0.0048154093097913355, 0.0048619350036057446, 0.0048154093097913355, 0.0045114840997070331, 0.0034906838696562147, 0.0040069911024866456, 0.0027766995669134334, 0.0016595801819374015, 0.0012182145074882836, 0.00098231827111984341, 0.00098231827111984363, 0.0012863691645616997, 0.0012395921040321833, 0.00093554121059032721, 0.0012629806342969417, 0.0010057068013846018, 0.0006081017868837127, 0.00032743942370661445, 4.6777060529516312e-05, 7.0165590794274467e-05, 7.0165590794274467e-05, 4.6777060529516745e-05]

## y_axis values must be normalised
sum_ys = sum(y_axis)

# normalize to 1
y_axis = [_/sum_ys for _ in y_axis]

# def gamma_f(x, a, loc, scale):
#     return gamma.pdf(x, a, loc, scale)

def norm_f(x, loc, scale):
#     print 'loc: ', loc, 'scale: ', scale, "\n"
    return norm.pdf(x, loc, scale)

fitting = norm_f

# param_bounds = ([-np.inf,0,-np.inf],[np.inf,2,np.inf])
result = curve_fit(fitting, x_axis, y_axis)
result_mod = result

# mod scale
# results_adj  = [result_mod[0][0]*.75, result_mod[0][1]*.85]

plt.plot(x_axis, y_axis, 'ro')
plt.bar(x_axis, y_axis, 1, alpha=0.75)
plt.plot(x_axis, [fitting(_, *result[0]) for _ in x_axis], 'b-')
plt.axis([0,35,0,.1])

# convert back into probability
y_norm_fit = [fitting(_, *result[0]) for _ in x_axis]
y_fit = [_*sum_ys for _ in y_norm_fit]
print list(y_fit)

plt.show()

我试图回答两个问题：

1. 这是我从正态分布曲线中得到的最佳拟合吗？我怎样才能提高我的合身度？

正态分布结果：

1. 如何为这些数据拟合对数正态分布，或者是否有更好的分布可以使用？

我在玩对数正态分布曲线调整 mu 和 sigma，看起来可能有更好的拟合。我不明白我做错了什么才能在 python 中得到类似的结果。

Answer 1

在 Python 中，我解释了一个技巧 here 如何使用 OpenTURNS 库非常简单地拟合 LogNormal：

import openturns as ot

n_times = [int(y_axis[i] * N) for i in range(len(y_axis))]
S = np.repeat(x_axis, n_times)

sample = ot.Sample([[p] for p in S])
fitdist = ot.LogNormalFactory().buildAsLogNormal(sample)

就是这样！

print(fitdist) 会告诉你>>> LogNormal(muLog = 2.92142, sigmaLog = 0.305, gamma = -6.24996)

而且看起来很合适：

import matplotlib.pyplot as plt

plt.hist(S, density =True, color = 'grey', bins = 34, alpha = 0.5)
plt.scatter(x_axis, y_axis, color= 'red')
plt.plot(x_axis, fitdist.computePDF(ot.Sample([[p] for p in x_axis])), color = 'black')
plt.show()

Answer 2

实际上，Gamma distribution 可能很适合@Glen_b 的建议。我正在使用带有 \alpha 和 \beta 的第二个定义。

注意：我用于快速拟合的技巧是计算均值和方差，对于典型的双参数分布，它足以恢复参数并快速了解它是否合适。

代码

import math
from scipy.misc import comb

import matplotlib.pyplot as plt

y_axis = [0.00032425299473065838, 0.00063714106162861229, 0.00027009331177605913, 0.00096672396877715144, 0.002388766809835889, 0.0042233337680543182, 0.0053072824980722137, 0.0061291327849408699, 0.0064555344006149871, 0.0065601228278316746, 0.0052574034010282218, 0.0057924488798939255, 0.0048154093097913355, 0.0048619350036057446, 0.0048154093097913355, 0.0045114840997070331, 0.0034906838696562147, 0.0040069911024866456, 0.0027766995669134334, 0.0016595801819374015, 0.0012182145074882836, 0.00098231827111984341, 0.00098231827111984363, 0.0012863691645616997, 0.0012395921040321833, 0.00093554121059032721, 0.0012629806342969417, 0.0010057068013846018, 0.0006081017868837127, 0.00032743942370661445, 4.6777060529516312e-05, 7.0165590794274467e-05, 7.0165590794274467e-05, 4.6777060529516745e-05]
x_axis = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0]

## y_axis values must be normalised
sum_ys = sum(y_axis)

# normalize to 1
y_axis = [_/sum_ys for _ in y_axis]

m = 0.0
for k in range(0, len(x_axis)):
    m += y_axis[k] * x_axis[k]

v = 0.0
for k in range(0, len(x_axis)):
    t = (x_axis[k] - m)
    v += y_axis[k] * t * t

print(m, v)

b = m/v
a = m * b

print(a, b)

z = []
for k in range(0, len(x_axis)):
    q = b**a * x_axis[k]**(a-1.0) * math.exp( - b*x_axis[k] ) / math.gamma(a)
    z.append(q)

plt.plot(x_axis, y_axis, 'ro')
plt.plot(x_axis, z, 'b*')
plt.axis([0, 35, 0, .1])
plt.show()

Answer 3

请注意，如果对数正态曲线是正确的并且您对两个变量都取对数，则应该有二次关系；即使这不是最终模型的合适尺度（因为方差效应——如果你的方差在原始尺度上接近恒定，它将超过小值），它至少应该为非线性拟合提供一个良好的起点。

确实，除了前两点，这看起来相当不错：

-- 对实心点的二次拟合可以很好地描述该数据，并且如果您想进行非线性拟合，应该给出合适的起始值。

（如果 x 中的错误完全可能，那么在最低 x 处的不拟合可能与 x 中的错误和 y 中的错误一样多）

顺便说一句，该图似乎暗示 gamma 曲线可能比对数正态曲线更适合整体（特别是如果您不想要减少前两点相对于第 4-6 点的影响）。通过在 x 和 log(x) 上回归 log(y) 可以得到一个很好的初始拟合：

缩放后的伽马密度为 g = cx^(a-1) exp(-bx) ... 取对数，得到 log(g) = log(c) + (a-1) log(x) - bx = b0 + b1 log(x) + b2 x ...因此将 log(x) 和 x 提供给线性回归例程将适合。关于方差效应的相同警告也适用（因此，如果您在 y 中的相对误差不是几乎恒定的，那么它可能最好作为非线性最小二乘拟合的起点）。

Answer 4

离散分布可能看起来更好 - 毕竟你的 x 都是整数。您的分布方差比平均值高约 3 倍，不对称 - 所以很可能像 Negative Binomial 这样的东西可能会很好地工作。这里是快速安装

r 略高于 6，因此您可能希望使用真正的 r 进行分发 - Polya 分发。

代码

from scipy.misc import comb

import matplotlib.pyplot as plt

y_axis = [0.00032425299473065838, 0.00063714106162861229, 0.00027009331177605913, 0.00096672396877715144, 0.002388766809835889, 0.0042233337680543182, 0.0053072824980722137, 0.0061291327849408699, 0.0064555344006149871, 0.0065601228278316746, 0.0052574034010282218, 0.0057924488798939255, 0.0048154093097913355, 0.0048619350036057446, 0.0048154093097913355, 0.0045114840997070331, 0.0034906838696562147, 0.0040069911024866456, 0.0027766995669134334, 0.0016595801819374015, 0.0012182145074882836, 0.00098231827111984341, 0.00098231827111984363, 0.0012863691645616997, 0.0012395921040321833, 0.00093554121059032721, 0.0012629806342969417, 0.0010057068013846018, 0.0006081017868837127, 0.00032743942370661445, 4.6777060529516312e-05, 7.0165590794274467e-05, 7.0165590794274467e-05, 4.6777060529516745e-05]
x_axis = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0]

## y_axis values must be normalised
sum_ys = sum(y_axis)

# normalize to 1
y_axis = [_/sum_ys for _ in y_axis]

s = 1.0 # shift by 1 to have them all at 0
m = 0.0
for k in range(0, len(x_axis)):
    m += y_axis[k] * (x_axis[k] - s)

v = 0.0
for k in range(0, len(x_axis)):
    t = (x_axis[k] - s - m)
    v += y_axis[k] * t * t

print(m, v)

p = 1.0 - m/v
r = int(m*(1.0 - p) / p)

print(p, r)

z = []
for k in range(0, len(x_axis)):
    q = comb(k + r - 1, k) * (1.0 - p)**r * p**k
    z.append(q)

plt.plot(x_axis, y_axis, 'ro')
plt.plot(x_axis, z, 'b*')
plt.axis([0, 35, 0, .1])
plt.show()