跳转至

多元回归和方差分析

多元线性回归

Python 中的线性代数

矩阵

import numpy as np
a = np.array([[1,2],[3,4]]) # 注意两个方括号

初始化矩阵

# 注意两个括号
A = np.zeros((2,2))
B = np.ones((2,2))
C = np.eye(2)   # 谐音 I
D = np.diag([1,2,3])
E = np.array(np.random.randint(10,size=(3,4)))  # 随机数矩阵

取值

A[1,1]  # a_22
A[1,:]   # 第二行
A[:,1]   # 第二列
A[:,1:3]    # 第二到第三列

矩阵运算

# 加法
A+B
# 每个元素的乘法(Element-wise)
A*B
# 每个元素的除法
A/B
  • Broadcast
# Broadcast
C = np.array(np.random.randint(3,size=(1,4)))
# C: 1x4, E: 3x4
C*E # C与E的每一行相乘

# F: 3x1
F = np.array(np.random.randint(3,size=(3,1)))

F*E # F与E的每一列相乘
矩阵转置
A.T
矩阵乘法
A = np.array(np.random.randint(10,size=(3,4)))
B = np.array(np.random.randint(10,size=(4,1)))

A@B
  • 查看矩阵的形状
A.shape
矩阵求逆
# 解 Ax = b
b = np.array(np.random.randint(10,size=(3,1)))
A = np.array(np.random.randint(10,size=(3,3)))

# 逆矩阵
A_inv = np.linalg.inv(A)
# 解方程
np.linalg.solve(A,b)
np.linalg.inv(A)@b  # 更慢

solve(A,b)inv(A)@b

solve更快,因为它使用了更高效的算法

计时比较:

%timeit np.linalg.solve(A,b)
%timeit np.linalg.inv(A)@b

多变量的回归

  • Y 由多个变量控制:
Y=b0+b1X1+b2X2+...+bpXp+ϵ

bkYxk 的回归系数.

  • 写为 β0,β1,...,βp 的形式:
Y=β0+β1X1+β2X2+...+βpXp+ϵ

多元线性函数的矩阵表示

  • 用列向量来表示 βY 的测量值
β=(β0β1βp),Y=(Y1Y2Yn)
  • 用矩阵表示 X 的测量值
X=(1x11x12x1p1xn1xn2xnp)

残差

δ=(y1β0x11β1x1pβpynβ0xn1β1xnpβp)=YXβ
δ2=(YXβ)T(YXβ)

最小二乘法求解

  • 求解 β 使得 δ2 最小
δ2β=0
  • 矩阵求导(x 为列向量)
xTax=aTxx=a
xTAxx=(A+AT)x

一顿推导TAT...

得到

XTXβ=XTY
β=(XTX)1XTY

L=(XTX),σ2L1β 的协方差矩阵,且有

  • var(β0)=σ2/n
  • var(βi)=σ2Lii1
  • Cov(βi,βj)=σ2Lij1

  • σ2 仍然可用残差来估计

δ2=δTδs=(YXβ)T(YXβ)
σ2=δ2np1

βj 的区间估计

(β^jβj)/(σ^Ljj1)t(np1)
from sklearn import linear_model
import numpy as np

# x,y 为 1*20 的行向量

X = np.vstack([x,y]).T      # 20*2
ols = linear_model.LinearRegression()   # 创建线性回归对象. OLS: Ordinary Least Squares
model = ols.fit(X, z)    # 拟合

x_pred = np.linspace(0, 10, 50)  # 生成预测值
y_pred = np.linspace(0, 8, 50)
x_predmesh, y_predmesh = np.meshgrid(x_pred, y_pred)  # 生成网格

# flatten() 将多维数组转换为一维数组
model_viz = np.array([x_predmesh.flatten(), y_predmesh.flatten()]).T    # 2500*2
print(model_viz.shape)
predicted = model.predict(model_viz)

fig = plt.figure(figsize=(4, 4))
ax = fig.add_subplot(111, projection='3d')
ax.plot(x_predmesh.flatten(), y_predmesh.flatten(), predicted, , alpha=0.9)

ax.scatter(x, y, z, color='red')

矩阵运算求线性回归系数

ones = np.ones(20)
Xm = np.vstack([ones, x, y]).T
Y = np.array(z).reshape(20, 1)
b = np.linalg.solve(Xm.T @ Xm, Xm.T @ Y)
print(b)

矩阵运算求线性回归误差

statsmodels 无脑求解:

import statsmodel.api as sm

Xp = sm.add_constant(X) # Xp = np.vstack([ones, x, y]).T

# OLS 拟合
model = sm.OLS(z, Xp)
Fit_result = model.fit()
print(Fit_result.summary())

多元线性回归的假设检验

  • 单个系数 βi 是否为0,由 βit 分布得到:

\((β^iβi)/(σ^Lii1)t(np1)\)

  • 一系列 β1,,βr(r<p) 为0,即 X1,,XrY 无显著影响
  • H0: β1==βr=0

\(R1=δ2\)

\(R2=δ2|β1==βr=0=i=1n(Yiβ0)2\)

  • δ2

转化为线性回归的模型

  • 多项式回归
Y=β0+β1x+β2x2++βpxp+ϵ
  • 指数拟合
from sklearn.preprocessing import PolynomialFeatures

a = 0.3
b = -4
c = 10
d = 8

x2 = st.uniform(0, 10).rvs(40)
x2 = np.sort(x2)
y2 = 

偏相关

  • 一组观测量 Xi 中,相互之间有关联
  • X3 为一个人的收入,X1,X2 为ta在吃和穿上的花费
  • 观察到 X1,X2 之间的正相关
  • 但实际上 X1,X2 均由 X3 所带动,使得它们呈现正相关
  • 若能去掉 X3 的影响,X1,X2 之间的实际相关性可能转为负相关
  • 重新定义 X1,X2

\(X1=X1L(X3)\)

\(X2=X2L(X3)\)

\(ρ(X1,X2)=ρ12ρ13ρ23(1ρ132)(1ρ232)\)

  • 将原本 X 的相关矩阵写为

\(ρ=[1ρ12ρ13ρ211ρ23ρ31ρ321]\)

Pij


大作业

自己找数据和主题!

参考数据来源: - Our World in Data


方差分析(ANOVA)

单因素

同因素的不同水平对结果的影响是否显著不同

  • 例:配方A,B,各自得到一组关于指标Y的数据
[Y11,Y12,,Y1n],[Y21,Y22,,Y2m]

配方A、B对结果的影响是否显著不同?

  • 计算Y在A、B配方下的平均值. 假设 H0:Y¯A=Y¯B
S2=1n+m2[i=1n(YAiY¯A)2+i=1m(YBiY¯B)2]
T=nmn+m(Y¯AY¯Bθ0)St(n+m2)
x=[2.74,2.75,2.72,2.69]
y=[2.75,2.78,2.74,2.76,2.72]
print(st.ttest_ind(x,y))

水平 k 多于 2 时

k 个水平,每个水平下有 ni 个数据,总数据量为 n=i=1kni.

水平 i 对 Y 的影响造成的平均值记为 ai,对于只受到 ai 影响的数据 Yij,有

Yij=ai+ϵij
  • 假设 H0:a1=a2==ak
  • 各水平间的差异: SSA=i=1kni(Y¯iY¯)2
  • 同水平内样本的差异: SSe=i=1kj=1ni(YijY¯i)2
  • 总样本差异: SS=i=1kj=1ni(YijY¯)2

可以证明

SS=SSA+SSe
F=SSA/(k1)SSe/(nk)F(k1,nk)
  • SSA/SSe<Fα(k1,nk) 时,假设 H0 成立
meanx = np.mean(x)
meany = np.mean(y)
meanall = np.float64(sum(x+y)/len(x+y))

ssa = len(x)*((meanx-meanall)**2)+len(y)*((meany-meanall)**2)
sse = sum((x-meanx)**2)+sum((y-meany)**2)
sta = ssa/(sse/(len(x)+len(y)-2))

print(sta,1-st.f(1,len(x)+len(y)-2).cdf(sta))
print(st.f_oneway(x,y)) # 单因素 F-test
print(st.ttest_ind(x,y))

可以发现三种方式得到的 p-value 是一样的!

多因素

当因素有多种时,分析每种因素对结果的影响

因素 A,B,各自有 k,l 个水平,将