0

下面我使用 nympy 生成的绘图来比较帕累托分布的经验矩和分析矩（我使用以下链接获取均值、方差、偏度和过度峰度的公式，https://en.wikipedia.org/wiki/Pareto_distribution）。获得的均值和方差结果相似，但偏度和峰度的结果却大不相同。我究竟做错了什么？

事先谢谢你

import numpy as np
import pandas as pd
from scipy.stats import skew
from scipy.stats import kurtosis
from prettytable import PrettyTable

x_m = [1, 2, 3, 4] 
alpha = [5, 6, 7, 8]

#drawing samples from distribution

for a, x in zip(alpha, x_m):
    print (a, x)
data = (np.random.default_rng().pareto(a, 10000000)+1) * x
mean = np.mean(data) 
var = np.var(data) 
skew = skew(data) 
kurt = kurtosis(data)

#Analytical estimation

for a, x in zip(alpha, x_m):
    a_mean = (a*x)/(a-1)
    a_var = (a*x**2)/((a-1)**2*(a-2))
    a_skew = (2*(1+a)/(a-3))*(np.sqrt(a-2/a))
    a_kurt = (6*(a**3+a**2-6*a-2))/(a*(a-3)*(a-4))
    
#Table

header = ['Moments', 'Simulated', 'Analytical']
Moments = ['Mean', 'Variance', 'Skewness', 'Excess Kurtosis']
Simulated = [round(mean,4), round(var,4), round(skew,4), round(kurt,4)]
Analytical = [round(a_mean,4), round(a_var,4), round(a_skew,4), round(a_kurt,4)]
table = PrettyTable()

table.add_column(header[0], Moments)
table.add_column(header[1], Simulated)
table.add_column(header[2], Analytical)
print(table)
4

1 回答 1

0

这是分析偏度线上的一个小错字：

a_skew = (2*(1+a)/(a-3))*(np.sqrt(a-2/a))

应该

a_skew = (2*(1+a)/(a-3))*(np.sqrt((a-2)/a))

如果样本足够大，峰度是正确的。

于 2020-11-15T13:18:34.047 回答