关于本文多元正态分布的推断理论,请参考《应用多元统计分析及带R的相关主题》第五章,其中提供了相关例子的R代码。对于偏爱Python的我来说,希望通过Python得到同样的结果。数据
关于本文多元正态分布的推断理论,请参考《应用多元统计分析及带R的相关主题》第五章,其中提供了相关例子的R代码。对于偏爱Python的我来说,希望通过Python得到同样的结果。
数据集的URL:
www.stat.ubc.ca/~·朗/text。
使用得数据集示例有:class.dat2,consum2000.txt,consume 2010 . txt .
Python编程分析前,先用R软件把数据集转换成格式。虽然Python也可以读取txt文件,但是我更喜欢读取csv格式,所以通过下面的代码,将数据集转换成CSV格式保存在本地。
data = read.table('class.dat2', header = T)write.csv(data, 'class.csv')
例1
需要类别数据集。这个例子可以简单概括为期中考试前有两次考试quiz1和quiz2,期中考试后有两次考试quiz3和quiz4,比较学生在quiz1和quiz2之间的成绩有没有进步,学生在quiz3和quiz4之间的成绩有没有进步。makeμ1 = mean(quiz 1 –quiz2),μ2 =均值(quiz3 & # 8211第四题.
推导多元正态分布,其编程思想是:
1.导入数据->: 2。求解样本均值和样本协方差矩阵->: 3。计算霍特林’的统计->: 4。计算P值,根据P值结合实例分析结果。
代码实现如下:
# 引入第三方库import pandas as pdimport numpy as npfrom scipy import stats# 读取数据data = pd.read_csv("class.csv")df = pd.DataFrame(data)# 构建矩阵y = np.c_[df.quiz1-df.quiz2, df.quiz3-df.quiz4]print(y)# 计算样本数量nn = np.shape(y)[0]# 计算变量数目pp = np.shape(y)[1]# 计算样本均值y = pd.DataFrame(y)y_bar = y.mean()print(y_bar)# 计算样本协方差S_y = y.cov()print(S_y)# 计算 Hotelling's T statisticT_sq = n * np.dot(np.dot(y_bar.T, np.linalg.inv(S_y)), y_bar)T_sq2 = ((n - p)/(p * (n - 1))) * T_sqprint('T_sq2:', T_sq2)# 计算p值p_value = 1 - stats.f.cdf(T_sq2, p, n-p)print('p_value:', p_value)
输出结果:
p_value: 0.05442091231270707
从p值可以看出,quiz1和quiz2之间以及quiz3和quiz4之间存在一些差异,但这些差异在5%的水平上不具有统计学意义。
例2
需要两个数据集,消费2000和消费2010。这个例子可以简单概括为2000年和2010年在食品、布料、居民(Resid)、交通(TranC)、教育(Educ)等消费结构上有没有变化的比较。制造
μ1 =平均值(Food.2010 & # 8211food . 2000);
μ2 =均值(Cloth.2010 & # 8211布. 2000);
μ3 =平均值(Resid.2010 & # 8211resid . 2000);
μ4 =平均值(TranC.2010 & # 8211tranc . 2000);
μ5 =平均值(Educ.2010 & # 8211Educ.2000)。
代码实现:
# 引入第三方库import numpy as npimport pandas as pdfrom scipy import stats# 导入数据consum00 = pd.read_csv("consum2000.csv")consum10 = pd.read_csv("consum2010.csv")# 计算2010年支出份额data10 = consum10.iloc[:, 1:9]sum10 = data10.sum(axis=1)X = data10.div(sum10, axis='rows')print(X)# 计算2000年支出份额data00 = consum00.iloc[:, 1:9]sum00 = data00.sum(axis=1)Y = data00.div(sum00, axis='rows')print(Y)# 求X与Y之差XY_d = np.c_[X.iloc[:, 0:3]-Y.iloc[:, 0:3], X.iloc[:, 5:7]-Y.iloc[:, 5:7]]XY_d = pd.DataFrame(XY_d, columns=('Food', 'Cloth', 'Resid', 'TranC', 'Educ'))# 计算样本均值d_mean = XY_d.mean()print(d_mean)# 计算样本协方差阵d_S = XY_d.cov()# 计算样本大小n = np.shape(XY_d)[0]# 计算变量数p = np.shape(XY_d)[1]# 计算 Hotelling's T 统计量T2 = n * np.dot(np.dot(d_mean.T, np.linalg.inv(d_S)), d_mean)Tstar2 = ((n-p)/(p*(n-1)))*T2# 计算p值p_value = 1 - stats.f.cdf(Tstar2, p, n-p)print('p_value:', p_value)
输出结果:
p_value: 7.460698725481052e-14
可以看出,P值近似为0,拒绝了原来的假设,说明2000年和2010年的消费结构发生了明显的变化。