如何使用 skewnorm 生成具有指定偏斜的分布？

Question

我正在尝试生成一个随机分布，在其中控制均值、SD、偏度和峰度。

生成分布后，我可以通过一些简单的数学计算来求解平均值和标准差。

峰度 我暂时将其搁置，因为它看起来太难了。

偏斜是今天的问题。

import scipy.stats

def convert_to_alpha(s):
    d=(np.pi/2*((abs(s)**(2/3))/(abs(s)**(2/3)+((4-np.pi)/2)**(2/3))))**0.5 
    a=((d)/((1-d**2)**.5))
    return(a)

for skewness_expected in (.5, .9, 1.3):
    alpha = convert_to_alpha(skewness_expected)
    r = stats.skewnorm.rvs(alpha,size=10000)
    print('Skewness expected:',skewness_expected)
    print('Skewness obtained:',stats.skew(r))
    print()

Skewness expected: 0.5
Skewness obtained: 0.47851348006629035

Skewness expected: 0.9
Skewness obtained: 0.8917020428586827

Skewness expected: 1.3
Skewness obtained: (1.2794406116842627+0.01780402125888404j)

我知道计算出的偏度通常与所需的偏度不匹配 - 毕竟这是一个随机分布。但我很困惑如何获得偏度 > 1 的分布而不陷入复数范围。rvs 方法似乎无法处理它，因为只要偏度 > 1，参数 alpha 就是一个虚数。

我该如何修复它，以便生成偏度 > 1 的分布，但不会出现复数？

[感谢Warren Weckesser为我指明了维基百科，以便我编写 Convert_to_alpha 函数。]

Answer 1

了解这个帖子已经有一年半的历史了，但我最近也遇到了这个问题，而且似乎从未在这里得到解答。在 stats.skewnorm 的 alpha 和偏度统计（顺便说一句，这是一个很好的函数）之间进行转换的另一个问题是，这样做也会改变分布的集中趋势的度量，这对我的需求来说是有问题的。

我是基于 F-distribution ( https://en.wikipedia.org/wiki/F-distribution ) 开发的。大量工作的最终结果是您指定所需的均值、SD 和偏度以及所需样本大小的函数。如果有人愿意，我可以分享其背后的工作。在极端设置下，输出标清和偏斜变得有点粗糙。大概是因为 F 分布自然位于 1 左右。当偏斜值接近于零时，这也是一个很大的问题，在这种情况下，无论如何都不需要这个函数。

from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

def createSkewDist(mean, sd, skew, size):

    # calculate the degrees of freedom 1 required to obtain the specific skewness statistic, derived from simulations
    loglog_slope=-2.211897875506251 
    loglog_intercept=1.002555437670879 
    df2=500
    df1 = 10**(loglog_slope*np.log10(abs(skew)) + loglog_intercept)

    # sample from F distribution
    fsample = np.sort(stats.f(df1, df2).rvs(size=size))

    # adjust the variance by scaling the distance from each point to the distribution mean by a constant, derived from simulations
    k1_slope = 0.5670830069364579
    k1_intercept = -0.09239985798819927
    k2_slope = 0.5823114978219056
    k2_intercept = -0.11748300123471256

    scaling_slope = abs(skew)*k1_slope + k1_intercept
    scaling_intercept = abs(skew)*k2_slope + k2_intercept

    scale_factor = (sd - scaling_intercept)/scaling_slope    
    new_dist = (fsample - np.mean(fsample))*scale_factor + fsample

    # flip the distribution if specified skew is negative
    if skew < 0:
        new_dist = np.mean(new_dist) - new_dist

    # adjust the distribution mean to the specified value
    final_dist = new_dist + (mean - np.mean(new_dist))

    return final_dist




'''EXAMPLE'''
desired_mean = 497.68
desired_skew = -1.75
desired_sd = 77.24

final_dist = createSkewDist(mean=desired_mean, sd=desired_sd, skew=desired_skew, size=1000000)

# inspect the plots & moments, try random sample
fig, ax = plt.subplots(figsize=(12,7))
sns.distplot(final_dist, hist=True, ax=ax, color='green', label='generated distribution')
sns.distplot(np.random.choice(final_dist, size=100), hist=True, ax=ax, color='red', hist_kws={'alpha':.2}, label='sample n=100')
ax.legend()

print('Input mean: ', desired_mean)
print('Result mean: ', np.mean(final_dist),'\n')

print('Input SD: ', desired_sd)
print('Result SD: ', np.std(final_dist),'\n')

print('Input skew: ', desired_skew)
print('Result skew: ', stats.skew(final_dist))

输入平均值：497.68
结果平均值：497.6799999999999

输入标清：77.24
结果标清：71.69030764848961

输入偏差：-1.75
结果偏差：-1.6724486459469905