python - 使用 SciPy 拟合 Levy-Stable 分布

Question

在 1.2 SciPy 中添加了 ability to fit a Levy-Stable distribution .我有一些我想要拟合的分布，但在运行拟合时遇到了一些问题。

这是我的测试用例:

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 Docs在这个问题上有点单薄。

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

Answer 1

Scipy 的 levy 稳定分布实现主要使用 Nolan 的方法，将参数空间(alpha，beta)分成几个部分，其中一些需要棘手的积分来评估。

Scipy 使用 MLE 估计参数，由于这些相同的积分，这可能非常慢。有用于评估 levy 稳定 PDF 的实验性 FFT 支持，此功能有望通过 this 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)

希望对您有所帮助。