问题描述
我不明白 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)