使用 SciPy 拟合具有正偏斜与负偏斜的 Levy-Stable 分布

如何解决使用 SciPy 拟合具有正偏斜与负偏斜的 Levy-Stable 分布

我不明白 scipy.stats.levy_stable 的 _fitstart() 方法返回的参数用于具有正和负 beta 参数的分布。直观地说，在生成随机样本时更改 beta 的符号不应影响拟合数据时对 alpha 的估计。我不确定 beta 的符号应该对 _fitstart() 返回的第三个参数有什么影响，但我希望在按照 this answer 的建议转换返回值后，符号可能会反转。

from scipy.stats import levy_stable
from scipy.stats import rv_continuous as rvc
import numpy as np

points = 1000000
jennys_constant = 8675309

pconv = lambda alpha,beta,mu,sigma: (alpha,mu - sigma * beta * np.tan(np.pi * alpha / 2.0),sigma)

rvc.random_state = jennys_constant

def test_fitstart(alpha,beta):
    draw = levy_stable.rvs(alpha,size=points)
    
    # use scipy's quantile estimator to estimate the parameters and convert to S parameterization
    return pconv(*levy_stable._fitstart(draw))

print("A few calls with beta=1")
for i in range(3):
    print(test_fitstart(alpha=1.3,beta=1))

print("A few calls with beta=-1")
for i in range(3):
    print(test_fitstart(alpha=1.3,beta=-1))

>>> A few calls with beta=1
>>> (1.3059810788754223,1.0,1.9212069030505312,1.0017497273563876)
>>> (1.3048494867305243,1.92631956349381,1.000064636906844)
>>> (1.3010492983811222,1.9544520781484407,0.9999042085058586)
>>> A few calls with beta=-1
>>> (1.3652389860952416,-1.0,0.3424825654388899,1.0317366391952136)
>>> (1.370069101697994,0.3560781956631771,1.0397745333221347)
>>> (1.3682310757082936,0.34621980810217745,1.037169706715312)

查看 _fitstart() 代码，我认为 alpha 的查找可能应该使用 nu_beta 的绝对值，但事实并非如此，因此查找可能是在 nu_beta_range 之外进行推断。

同样，我想知道是否应该在 delta 的计算中使用某些东西的绝对值，在应用裁剪之前，对 beta 的符号进行裁剪后调整？实际上，再看一遍我认为应该将裁剪应用于 c（缩放参数，必须为正）。不应将裁剪应用于 delta（位置参数 = 平均值，可从 -inf 到 inf）。对吗？

解决方法

levy_stable._fitstart() 没有正确处理负偏斜数据，但我们可以通过反映关于原点的样本来解决这个问题。 _fitstart() 然后将返回对稳定性和尺度参数的合理估计，这些参数不受反射影响。偏度和loc参数的估计在反射样本中被反演。

一个简单的包装函数可以在调用 _fitstart() 之前检查数据是向右还是向左倾斜，然后根据需要反转反转参数估计。这不会修复 levy_stable.fit() 本身，但至少我们可以从 _fitstart() 获得分位数估计。

import numpy as np
from scipy import __version__ as scipy_version
from scipy.stats import levy_stable

points = 1000000

const = 314

def lsfitstart(data):
    """Wrapper for levy_stable._fitstart() to fix data with negative skew"""
    skewleft = np.mean(data) <= np.median(data)
    
    # reverse sign of the data points if distribution has negative skew
    alpha,beta,loc,scale = levy_stable._fitstart(-data if skewleft else data)
    
    # reverse sign of skewness and loc estimates if distribution has negative skew
    beta_fixed,loc_fixed = [-x if skewleft else x for x in (beta,loc)]
    
    # clip scale parameter to ensure it is positive
    scale_fixed = np.clip(scale,np.finfo(float).eps,np.inf)

    return (alpha,beta_fixed,loc_fixed,scale_fixed)


print(scipy_version)

sample = levy_stable.rvs(alpha=1.3,beta=1,size=points,random_state=const)

print("levy_stable fit        : alpha (stabililty),beta (skewness),scale")
print("_fitstart(positive    ): ",levy_stable._fitstart(sample))
print("_fitstart(negative=bad): ",levy_stable._fitstart(-sample))
print()

print("lsfitstart(positive   ): ",lsfitstart(sample))
print("lsfitstart(negative=OK): ",lsfitstart(-sample))
>>> 1.5.2
>>> levy_stable fit        : alpha (stabililty),scale
>>> _fitstart(positive    ):  (1.3055214922752139,1.0,2.220446049250313e-16,1.0002159057207403)
>>> _fitstart(negative=bad):  (1.3673555827622812,-1.0,1.9389262857717497,1.0337386531320203)
>>> 
>>> lsfitstart(positive   ):  (1.3055214922752139,1.0002159057207403)
>>> lsfitstart(negative=OK):  (1.3055214922752139,-2.220446049250313e-16,1.0002159057207403)