使用 scipy 从 pdf 样本坐标创建最佳拟合概率分布答案

【问题标题】：Creating a best fit probability distribution from pdf sample coordinates with scipy使用 scipy 从 pdf 样本坐标创建最佳拟合概率分布
【发布时间】：2017-09-04 12:05:34
【问题描述】：

问题：我有数据点表示从概率分布采样的坐标（在这种情况下，我们将假设一个离散的概率分布函数）我们基本上是从 pdf 形成“pdf 的最佳拟合”数据在这里。

给定： pdf 的样本坐标和适合它的 pdf 类型的类型（例如 lognorm）

返回：理想情况下是 pdf 参数，或者是最佳拟合分布的坐标。

我没有在 stackoverflow 上找到这个问题/答案的问题，我知道这可能是不好的做法。似乎 scipy 明确喜欢从原始数据构建 pdf 参数，而不是来自 pdf 的样本坐标。

我有向量：

x = list(range(40))

y = 
[0.032935611986072325,
 0.15399668951796566,
 0.19217568076280733,
 0.16189644686218774,
 0.11504756998080325,
 0.09474568682103104,
 0.08971162676825704,
 0.06198299715985481,
 0.04408241680044377,
 0.026817519111333753,
 0.013562814925870696,
 0.007007365243147507,
 0.003909173588759217,
 0.0015053452905258473,
 0.00037481359597322736,
 0.0001378624720821066,
 5.734365756863486e-05,
 2.9711739672867803e-05,
 8.022169711674307e-06,
 5.942347934573561e-06,
 2.228380475465085e-06,
 3.7139674591084754e-06,
 8.913521901860341e-07,
 8.913521901860341e-07,
 5.94234793457356e-07,
 2.97117396728678e-07,
 2.97117396728678e-07,
 2.97117396728678e-07,
 1.48558698364339e-07,
 0.0,
 0.0,
 0.0,
 0.0,
 0.0,
 0.0,
 0.0,
 0.0,
 0.0,
 0.0]

【问题讨论】：

标签： python scipy

【解决方案1】：

致电您的 PDF f(x):

如果您的数据确实代表{x, f(x)}，那么您可以尝试简单地优化f 的参数，例如使用https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html#scipy.optimize.leastsq

另一方面，如果您的数据是来自概率分布的样本，即您的数据看起来像{x}，但每个x 的选择概率为f(x)，那么您应该尝试马尔科夫链蒙特卡罗来估计@987654330 @。 Python 有多种选择：

https://pystan.readthedocs.io/en/latest/

http://docs.pymc.io/notebooks/getting_started.html#Model-fitting

【讨论】：

【解决方案2】：

我认为您的数据代表自 sum(y) = 1 以来的 pdf {x, y = pdf(x)}。当我们绘制您的数据并稍加修正x = list(range(39)) 时，我们会得到一条类似于对数正态 (?) 的曲线。

import matplotlib.pyplot as plt

x = list(range(39))
plt.plot(x, y)

您可以用来避免优化算法的一个技巧是将数据转换为样本，因为每个 y[i] 与 x[i] 的频率成正比。换句话说，如果你想要一个大小为 N 的“完美”样本S，每个x[i] 将出现N * y[i] 次。

N = 20.000
n_times = [int(y[i] * N) for i in range(len(y))]
S = np.repeat(x, n_times)

剩下要做的就是将 LogNormal 分布拟合到 S。就我个人而言，我习惯于 OpenTURNS 库。您只需将S 格式化为ot.Sample，只需将其整形为维度为 1 的 N 个点

import openturns as ot

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

fitdist是一个“ot.Distribution”，你可以打印看看它的参数

print(fitdist)
>>> LogNormal(muLog = 1.62208, sigmaLog = 0.45679, gamma = -1.79583)

或使用fitdist.computePDF 以ot.Sample 格式作为参数的内置方法绘制两条曲线

plt.plot(x, y)
plt.plot(x, fitdist.computePDF(ot.Sample([[p] for p in x])))

【讨论】：