使用SciPy拟合Levy-Stable分布

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文档在这个主题上有点薄弱。

我遇到了类似的问题，难以拟合我的真实数据。

- rhaskett

---

也许你正在使用 Jenny 常数的值，链接。看起来你可能错过了一个小数点。希望这可以帮到你。 - TavoGLC

@TavoGLC，使用Jenny的常数作为random_state似乎有些幽默。但也许你的评论也是这样的。 :) - Ulrich Stern

Answer 1

Scipy的Levy稳定分布实现主要使用Nolan的方法将参数空间（alpha，beta）分成几个部分，其中一些需要评估棘手的积分。Scipy使用MLE估计参数，由于这些相同的积分，这可能非常慢。有关评估Levy稳定PDF的实验性FFT支持，希望此功能在1.3里程碑中用this PR完成后得到显着改进。但是，即使使用FFT，fit()方法仍然相当缓慢。有一个更快的分位数估计器（McCulloch），用作分布参数的第一次猜测（在使用fit()进行估计时）。可以直接使用_fitstart()调用它。话虽如此，似乎用于生成Scipy随机样本（从rvs()）的参数化与用于生成pdf / cdf的参数化不同。我希望将来能够研究一下。在那之前（如@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)

希望有所帮助。

- Blair Azzopardi

---

是的，所有其他的分布拟合得相当快，但是levy_stable要花费数小时的时间。 - NoName

@NoName 我认为缓慢是由于scipy的MLE方法拟合分布所导致的。在https://github.com/scipy/scipy/pull/9523中有一些优化和修复等待着，但没有人真正能够（或有时间）审查PR。 - Blair Azzopardi

Answer 2

看起来您已经正确设置了问题; 文档对于levy_stable的超类rv_continuous，包含了所有函数的链接（例如fit()）。 我的直觉是运行时间非常慢是SciPy的一个错误。

使用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)

- Ulrich Stern

---

有趣，pylevy在conda里吗？ - rhaskett

我的直觉是不是这样的；我不使用conda。我按照我的答案中链接的pylevy的GitHub页面上的安装说明进行操作。该页面有版本1.0，而PyPI上的版本相当旧（0.3）。 - Ulrich Stern