【Rで最小2乗フィット】QR分解 (Householder変換)を使った最小2乗法

北川先生の本「Rによる時系列モデリング入門」で，Householder変換を使った最小2乗フィットの話が登場したので，Rでやってみます．

逆行列を使った最小2乗法
QR分解 (Householder変換)を使った最小2乗法
まとめ
【付録】3変数の場合

　ここでは，以下の式を使ってデータを生成します．

　 $\displaystyle y_n = \sum_{i=1}^{2} a_i x_{ni} + \varepsilon_n$

ここで， $\displaystyle \varepsilon_n$ は，平均0，標準偏差sd.epsの正規乱数とします．

　以下のスクリプトでは，データの長さをN，パラメタ $a_1$ をa1， $a_2$ をa2に設定します．

　 $\vec{y} = \begin{bmatrix} y_1 \\ y_{2} \\ \vdots\\ y_{N} \\ \end{bmatrix}, \quad Z = \begin{bmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \\ \vdots & \vdots \\ x_{N1} & x_{N2} \\ \end{bmatrix}, \quad \vec{a} = \begin{bmatrix} a_1 \\ a_2 \\ \end{bmatrix} , \quad \vec{\varepsilon} = \begin{bmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \vdots \\ \varepsilon_{N} \\ \end{bmatrix}$

と表せば，

　 $\vec{y} = Z \vec{a} + \vec{\varepsilon}$

です．ベクトルはすべて縦書きの列ベクトルとします．

逆行列を使った最小2乗法

　残差としての $\{\varepsilon_{n}\}$ の2乗和が最小になるように $\vec{a }$ を計算すると，

　 $\vec{a } = \left({}^t Z Z \right)^{-1} {}^t Z \, \vec{y}$

となります． ${}^t Z$ は， $Z$ の転置行列を表します．説明は，過去の記事を参考にしてください．
chaos-kiyono.hatenablog.com

　以下は，N = 10，a1 = 2.4，a2 = -3.0とした例です．Rのコマンドを使って逆行列を計算して，パラメタを推定しています．

# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
##############################
# a.hat = (t(Z) Z)^{-1} t(Z) Y 
##############################
a.hat <- solve(t(Z) %*% Z) %*% (t(Z) %*% Y)
# 残差の２乗平均
sig2.hat <- sum((Y - Z %*% a.hat)^2)/nrow(Y)
# a1, a2の推定結果
a.hat[,1]
# 残差の2乗平均 sig2
sig2.hat

　このスクリプトを実行すると，

> # a1, a2の推定結果
> a.hat[,1]
[1]  2.425303 -3.634048
> # 残差の2乗平均 sig2
> sig2.hat
[1] 1.607804

のように，パラメタと残差が計算されます．

QR分解 (Householder変換)を使った最小2乗法

　次は，QR分解 (Householder変換)でパラメタを推定してみます．pracmaパッケージを使いました．最初に，このパッケージをインストールしてください．

　使ったRスクリプトは以下です．

#install.packages("pracma")
require(pracma)
# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
###
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# Householder変換
UX <- householder(X)$R
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat

　このスクリプトを実行すると，

> # a1, a2の推定結果
> a.hat
[1]  2.611575 -2.950162
> # 残差の2乗平均 sig2
> sig2.hat
[1] 0.6768868

のように，パラメタと残差が計算されます．

　pracmaパッケージを使わなくても，Rのコマンドqrを使えば，QR分解できます．中身のアルゴリズムについては，何を使っているのか分かりません．Rスクリプトは以下です．

# データの長さ
N <- 10
# 誤差の標準偏差
sd.eps <- 1
eps <- rnorm(N,0,sd.eps)
##############################
a1 <- 2.4
a2 <- -3.0
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N]),　ncol=2)
###
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# QR分解でRを求める
UX <- qr.R(qr(X))
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat

まとめ

　QR分解について，勉強してみてください．QR分解の方法には，グラム・シュミットの正規直交化法とか，Householder変換を使う方法とか，ギブンス回転を使う方法などがあるそうです．

【付録】3変数の場合

　逆行列を使う場合．

N <- 10
sig2 <- 1
eps <- rnorm(N,0,sqrt(sig2))
##############################
a1 <- 2.4
a2 <- -3.0
a3 <- -1.8
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + a3 * x3[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
x3 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+a3*x3+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N],x3[1:N]),　ncol=3)
##############################
# a.hat = (t(Z) Z)^{-1} t(Z) Y 
##############################
a.hat <- solve(t(Z) %*% Z) %*% (t(Z) %*% Y)
# 残差の２乗平均
sig2.hat <- sum((Y - Z %*% a.hat)^2)/nrow(Y)
# a1, a2の推定結果
a.hat[,1]
# 残差の2乗平均 sig2
sig2.hat

　Householder変換を使う場合．

#install.packages("pracma")
require(pracma)
###############
N <- 10
sig2 <- 1
eps <- rnorm(N,0,sqrt(sig2))
##############################
a1 <- 2.4
a2 <- -3.0
a3 <- -1.8
##############################
# 真のモデルは
# y[i] = a1 * x1[i] + a2 * x2[i] + a3 * x3[i] + eps[i]
##############################
x1 <- rnorm(N,0,1)
x2 <- rnorm(N,0,1)
x3 <- rnorm(N,0,1)
y <- a1*x1+a2*x2+a3*x3+eps 
##############################
Y <- matrix(y[1:N],　ncol=1) # ncolで列数のみ指定
Z <- matrix(c(x1[1:N],x2[1:N],x3[1:N]),　ncol=3)
##############################
# X = [Z|Y] と結合
X <- cbind(Z,Y)
##############################
# Householder変換
UX <- householder(X)$R
###
m <- ncol(UX)-1
a.hat <- numeric(m)
a.hat[m] <- UX[m,m+1]/UX[m,m]
for(i in 2:m){
 a.hat[m+1-i] <- (UX[m+1-i,m+1]-a.hat[(m+2-i):m] %*% UX[m+1-i,((m+2-i)):m])/UX[m+1-i,m+1-i]
}
# 残差の２乗平均
sig2.hat <- UX[3,3]^2/nrow(Y)
# a1, a2の推定結果
a.hat
# 残差の2乗平均 sig2
sig2.hat