多元正态分布的定义（多元统计分析经典例题）

关于本文多元正态分布推断（inference for a multivariate normal population）理论参见《applied multivariate statistical analysis and related topics with r》第五章内容，书中提供了相关示例的r代码，对于偏爱python的我，希望通过python得到同样的结果。

数据集的获取网址：
www.stat.ubc.ca/~lang/text.

示例用到的数据集分别为：class.dat2, consum2000.txt, consum2010.txt.

进行python编程分析前，先把数据集通过r软件转换下格式，虽然python也可以读取txt文件，但我更喜欢读取csv格式，所以通过以下代码，将数据集转换为csv格式并保存本地。

data = read.table('class.dat2', header = t)
write.csv(data, 'class.csv')

示例1

需要用到class数据集，这个示例可以简单概括为期中考试前有两次测试quiz1和quiz2，期中考试后有两次测试quiz3和quiz4，比较quiz1和quiz2之间学生成绩有无进步，以及quiz3和quiz4之间学生成绩有无进步。令μ1 = mean(quiz1 – quiz2)， μ2 = mean(quiz3 – quiz4)。

进行多元正态分布推断，编程思路为：

1.导入数据 -> 2.求解样本均值和样本协方差阵 -> 3.计算hotelling’s t 统计量 -> 4.计算p值，根据p值结合实例分析结果。

代码实现如下：

# 引入第三方库
import pandas as pd
import numpy as np
from 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)
# 计算样本数量n
n = np.shape(y)[0]
# 计算变量数目p
p = 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 statistic
t_sq = n * np.dot(np.dot(y_bar.t, np.linalg.inv(s_y)), y_bar)
t_sq2 = ((n - p)/(p * (n - 1))) * t_sq
print('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

需要用到consum2000, consum2010两个数据集，这个示例可以简单概括为比较2000年和2010年在食品（food）、衣物（cloth）、居民数（resid）、交通（tranc）以及教育（educ）消费结构有无变化。令

μ1 = mean(food.2010 – food.2000);

μ2 = mean(cloth.2010 – cloth.2000);

μ3 = mean(resid.2010 – resid.2000);

μ4 = mean(tranc.2010 – tranc.2000);

μ5 = mean(educ.2010 – educ.2000).

代码实现：

# 引入第三方库
import numpy as np
import pandas as pd
from 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年的消费结构发生了明显的变化。