使用 SciPy 拟合征费稳定分布

Question

在 1.2 SciPy 添加了拟合 Levy-Stable 分布的能力。我有一些我想适合的发行版，但我在让适合运行时遇到了一些问题。

这是我的测试用例：

points = 1000
jennys_constant = 8675309
alpha, beta = 1.8, -0.5

draw = levy_stable.rvs(alpha, beta, size=points, random_state=jennys_constant)
print(levy_stable.fit(draw))

我觉得如果我从 Levy-Stable 分布中抽取，我应该能够很容易地适应该抽取。但是，我收到了很多如下警告，并且问题在 1000 点上花费了很长时间。

C:\anaconda3\lib\site-packages\scipy\stats\_continuous_distns.py:3857: IntegrationWarning: The integral is probably divergent, or slowly convergent.
intg = integrate.quad(f, -xi, np.pi/2, **intg_kwargs)[0]

我是否错误地设置了问题？的SciPy的文档是关于该主题的位薄。

我在拟合我的真实世界数据时遇到了类似的问题。

Answer 1

Scipy 对 levy 稳定分布的实现主要使用 Nolan 的方法，该方法将参数空间（alpha、beta）分成几个部分，其中一些部分需要计算复杂的积分。

Scipy 使用 MLE 估计参数，由于这些相同的积分，这可能会非常慢。有用于评估征费稳定 PDF 的实验性 FFT 支持，此功能有望随着此 PR的 1.3 里程碑而显着改进。但是，即使使用 FFT，fit() 方法似乎仍然很慢。

有一个更快的分位数估计器 (McCulloch) 用作分布参数的第一次猜测（使用 fit() 进行估计时）。这可以使用 _fitstart() 直接调用。

也就是说，用于生成 Scipy 随机样本（来自 rvs()）的参数化似乎与用于生成 pdfs/cdfs 的参数化不同。我希望将来能看到的东西。

在那之前（正如@Ulrich 在他们的回答中所建议的那样）您可以使用 pylevy 或使用 _fitstart() 来估计参数并在之后转换参数化。

from scipy.stats import levy_stable
import numpy as np

points = 1000000
jennys_constant = 8675309
alpha, beta = 1.8, -0.5

draw = levy_stable.rvs(alpha, beta, size=points, random_state=jennys_constant)

# use scipy's quantile estimator to estimate the parameters and convert to S parameterization
pconv = lambda alpha, beta, mu, sigma: (alpha, beta, mu - sigma * beta * np.tan(np.pi * alpha / 2.0), sigma)
pconv(*levy_stable._fitstart(draw))

>>> (1.7990380668349146, -0.5661063359664303,
      -0.012873575589969821, 0.998276003705684)

希望有帮助。

Answer 2

看来您已经正确设置了问题； 的超类 的文档具有其所有函数的链接（例如，）。我的预感是，运行速度非常慢是 SciPy 的 bug。rv_continuouslevy_stablefit()

使用pylevy似乎fit_levy()有效：

import scipy.stats as st, levy

points = 1000
jennys_constant = 8675309
alpha, beta = 1.8, -0.5

draw = st.levy_stable.rvs(alpha, beta, size=points, random_state=jennys_constant)
print(levy.fit_levy(draw))

结果看起来相当不错（而且fit_levy()速度相当快）：

(par=0, alpha=1.84, beta=-0.29, mu=0.11, sigma=1.00, 1863.61502664704)