PH525x series - Linear Algebra E
本章节将通过举例说明在进行数据分析时矩阵代数所起到的作用，以及如何利用矩阵代数方程构建线性模型、计算最小二乘估计。
在开始正文之前，先来补充两个概念：
-
单位矩阵（identity matrix）
形如的矩阵叫做单位矩阵：
单位矩阵与数字1相似：如果你用另一个矩阵与单位矩阵相乘，那么所得结果将与原矩阵相同。在R中，可使用diag(n)函数生成n*n的单位矩阵。 -
逆矩阵（inverse matrix）
设是一个n阶矩阵，那么
的逆矩阵记为
，并且
。并不是所有的矩阵都有逆矩阵，比如说，一个2行2列全部是1的矩阵就没有逆矩阵。
平均值
我们都知道样本均值的计算公式是，那么矩阵代数是如何计算平均值的呢？
首先，在矩阵代数中，要计算均值，首先要定义一个数值全是1的列的矩阵：
其次，根据矩阵相乘法则便可以推导出均值的计算过程：
在R中举例说明如下：
library(UsingR)
y <- father.son$sheight
print(mean(y))
## [1] 68.68407
N <- length(y)
Y<- matrix(y,N,1)
A <- matrix(1,N,1)
barY=t(A)%*%Y / N
print(barY)
## [,1]
## [1,] 68.68407
在R中还可以使用crossprod()函数计算矩阵与转置矩阵之间的乘积：
barY=crossprod(A,Y) / N
print(barY)
## [,1]
## [1,] 68.68407
方差
假设 :
那么方差便等于:
在R中，如果crossprod()中仅输入一个矩阵，默认计算该矩阵与其转置矩阵的乘积：
r <- y - barY
crossprod(r)/N
## [,1]
## [1,] 7.915196
library(rafalib)
popvar(y)
## [1] 7.915196
可以看到使用上述两种方法计算出的方差是相等的
线性模型
仍拿父与子身高问题举例说明，不清楚背景的可以去这里了解下。
假设：
那么模型便可以被写成：
更简单的一种写法是：
根据最小二乘法方程我们知道，求，其实就是求
（RSS）的最小值，也就是：
的最小值。
利用微积分求RSS的最小值
在矩阵代数中，有很多法则可供我们计算偏导方程，将导数设为0然后对求解，就可以解决最小二乘估计问题，上述方程的导数是
（
这块儿没有搞明白），若设其值等于0：
经过变换可得：
最终：
这样对未知参数的求解就变成了已知的自变量与因变量之间的计算，注意最小二程方程与的类似，其导数是
。
在R中计算LSE
library(UsingR)
x=father.son$fheight
y=father.son$sheight
X <- cbind(1,x)
betahat <- solve( t(X) %*% X ) %*% t(X) %*% y
###or
betahat <- solve( crossprod(X) ) %*% crossprod( X, y )
print(betahat)
# [,1]
# 33.886604
#x 0.514093
求出后，我们便可以带入任意
求
了：
newx <- seq(min(x),max(x),len=100)
X <- cbind(1,newx)
fitted <- X%*%betahat
plot(x,y,xlab="Father's height",ylab="Son's height")
lines(newx,fitted,col=2)
fitted.png
是在数据分析中最常用的几个公式之一，它可以应用在许多场景中，比如之前提到过的自由落体问题：
set.seed(1)
g <- 9.8 #meters per second
n <- 25
tt <- seq(0,3.4,len=n) #time in secs, t is a base function
d <- 56.67 - 0.5*g*tt^2 + rnorm(n,sd=1)
X <- cbind(1,tt,tt^2)
y <- d
betahat <- solve(crossprod(X))%*%crossprod(X,y)
print(betahat)
# [,1]
# 56.5317368
#tt 0.5013565
# -5.0386455
newtt <- seq(min(tt),max(tt),len=100)
X <- cbind(1,newtt,newtt^2)
fitted <- X%*%betahat
plot(tt,y,xlab="Time",ylab="Height")
lines(newtt,fitted,col=2)
falling.fitted.png
lm()函数
上述利用矩阵代数求解的过程可以很简单的用lm()函数代替，使用方法如下：
X <- cbind(tt,tt^2)
fit=lm(y~X)
summary(fit)
##
## Call:
## lm(formula = y ~ X)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5295 -0.4882 0.2537 0.6560 1.5455
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.5317 0.5451 103.701 <2e-16 ***
## Xtt 0.5014 0.7426 0.675 0.507
## X -5.0386 0.2110 -23.884 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9822 on 22 degrees of freedom
## Multiple R-squared: 0.9973, Adjusted R-squared: 0.997
## F-statistic: 4025 on 2 and 22 DF, p-value: < 2.2e-16
从上面可以看出，所得结果与我们之前计算是一致的。
