基于python进行抽样分布描述及实践详解

2019-11-25 11:47:22

字体：大中小

来源：转载

供稿：网友

本次选取泰坦尼克号的数据，利用python进行抽样分布描述及实践。

备注：数据集的原始数据是泰坦尼克号的数据，本次截取了其中的一部分数据进行学习。Age:年龄，指登船者的年龄。Fare:价格，指船票价格。Embark:登船的港口。

1、按照港口分类，使用python求出各类港口数据年龄、车票价格的统计量（均值、方差、标准差、变异系数等）。

import pandas as pddf = pd.read_excel('/Users/Downloads/data.xlsx',usecols = [1,2,3] )#拿到港口'Embarked'、年龄'Age'、价格'Fare'的数据df2 = df.groupby(['Embarked'])#按照港口'Embarked'分类后，查看 年龄、车票价格的统计量。# 变异系数 = 标准差/平均值def cv(data):  return data.std()/data.var()df2 = df.groupby(['Embarked']).agg(['count','min','max','median','mean','var','std',cv])df2 = df2.apply(lambda x:round(x,2))df2_age = df2['Age']df2_fare = df2['Fare']

分类后年龄及价格统计量描述数据如下图：

年龄统计量

价格统计量

2、画出价格的分布图像，验证数据服从何种分布（正态？卡方？还是T?）

2.1 画出船票的直方图：

plt.hist(df['Fare'],20,normed=1, alpha=0.75)plt.title('Fare')plt.grid(True)

船票价格的直方图及概率分布

2.2 验证是否符合正态分布？

#分别用kstest、shapiro、normaltest来验证分布系数ks_test = kstest(df['Fare'], 'norm')#KstestResult(statistic=0.99013849978633, pvalue=0.0)shapiro_test = shapiro(df['Fare'])#shapiroResult(0.5256513357162476, 7.001769945799311e-40)normaltest_test = normaltest(df['Fare'],axis=0) #NormaltestResult(statistic=715.0752414548335, pvalue=5.289130045259168e-156)

以上三种检测结果表明 p<5%,因此船票数据不符合正态分布。

绘制拟合正态分布曲线:

fare = df['Fare']plt.figure()fare.plot(kind = 'kde')   #原始数据的正态分布M_S = stats.norm.fit(fare)  #正态分布拟合的平均值loc，标准差 scalenormalDistribution = stats.norm(M_S[0], M_S[1])  # 绘制拟合的正态分布图x = np.linspace(normalDistribution.ppf(0.01), normalDistribution.ppf(0.99), 100)plt.plot(x, normalDistribution.pdf(x), c='orange')plt.xlabel('Fare about Titanic')plt.title('Titanic[Fare] on NormalDistribution', size=20)plt.legend(['Origin', 'NormDistribution'])

船票拟合正态分布曲线

2.3 验证是否符合T分布？

T_S = stats.t.fit(fare)df = T_S[0] loc = T_S[1] scale = T_S[2] x2 = stats.t.rvs(df=df, loc=loc, scale=scale, size=len(fare))D, p = stats.ks_2samp(fare, x2) # (0.25842696629213485 2.6844476044528504e-21)

p = 2.6844476044528504e-21 ，p < alpha，拒绝原假设，价格数据不符合t分布。

对票价数据进行T分布拟合：

plt.figure()fare.plot(kind = 'kde') TDistribution = stats.t(T_S[0], T_S[1],T_S[2])  # 绘制拟合的T分布图x = np.linspace(TDistribution.ppf(0.01), TDistribution.ppf(0.99), 100)plt.plot(x, TDistribution.pdf(x), c='orange')plt.xlabel('Fare about Titanic')plt.title('Titanic[Fare] on TDistribution', size=20)plt.legend(['Origin', 'TDistribution'])

票价拟合T分布

2.4 验证是否符合卡方分布？

chi_S = stats.chi2.fit(fare)df_chi = chi_S[0] loc_chi = chi_S[1] scale_chi = chi_S[2] x2 = stats.chi2.rvs(df=df_chi, loc=loc_chi, scale=scale_chi, size=len(fare))Df, pf = stats.ks_2samp(fare, x2) # (0.16292134831460675, 1.154755913291936e-08)

p = 1.154755913291936e-08 ，p < alpha，拒绝原假设，价格数据不符合卡方分布。

对票价数据进行卡方分布拟合

plt.figure()fare.plot(kind = 'kde') chiDistribution = stats.chi2(chi_S[0], chi_S[1],chi_S[2])  # 绘制拟合的正态分布图x = np.linspace(chiDistribution.ppf(0.01), chiDistribution.ppf(0.99), 100)plt.plot(x, chiDistribution.pdf(x), c='orange')plt.xlabel('Fare about Titanic')plt.title('Titanic[Fare] on chi-square_Distribution', size=20)plt.legend(['Origin', 'chi-square_Distribution'])

票价拟合卡方分布

3、按照港口分类，验证S与Q两个港口间的价格之差是否服从某种分布

S_fare = df[df['Embarked'] =='S']['Fare']Q_fare = df[df['Embarked'] =='Q']['Fare']C_fare = df[df['Embarked'] =='C']['Fare']S_fare.describe()count  554.000000mean   27.476284std    36.546362min    0.00000025%    8.05000050%    13.00000075%    27.862500max   263.000000Q_fare.describe()count  28.000000mean   18.265775std   21.843582min    6.75000025%    7.75000050%    7.75000075%   18.906250max   90.000000C_fare.describe()count  130.000000mean   68.296767std    90.557822min    4.01250025%    14.45420050%    36.25210075%    81.428100max   512.329200

按照港口分类后，S港口样本数<=554,Q港口样本数<=28,C港口样本数<=130。

总体不服从正态分布，所以需要当n比较大时，一般要求n>=30，两个样本均值之差的抽样分布可近似为正态分布。X2的总体容量为28，其样本容量不可能超过30，故其S港和Q港两个样本均值之差（E(X1)-E(X2)）的抽样分布不服从正态分布。

S港和C港两个样本均值之差（E(X1)-E(X3)）的抽样分布近似服从正态分布，其均值和方差分别为E(E(X1) - E(X3)) = E(E(X1)) - E(E(X3)) = μ1 - μ3；D(E(X1) + E(X3)) = D(E(X1)) + D(E(X3)) = σ1²/n1 + σ3²/n3 。绘图如下：

miu = np.mean(S_fare) - np.mean(C_fare)sig = np.sqrt(np.var(S_fare, ddof=1)/len(S_fare) + np.var(C_fare, ddof=1)/len(C_fare))x = np.arange(- 110, 50)y = stats.norm.pdf(x, miu, sig)plt.plot(x, y)plt.xlabel("S_Fare - C_Fare")plt.ylabel("Density")plt.title('Fare difference between S and C')plt.show()