多元統計分析經典例題：多元正態分布的定義

關於本文多元正態分布推斷（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)

輸出結果：