演習問題 2 の解答

データサイエンス Lec05 ^{演習問題解答} 演習問題 1 ^の解答

定義より

Cov[ ˆw] =E[( ˆw−E[ ˆw])( ˆw−E[ ˆw])^⊤].

ここで，

wˆ−E[ ˆw] = (X^⊤X)⁻¹X^⊤y−w

= (X^⊤X)⁻¹X^⊤(Xw+ε)−w

= (X^⊤X)⁻¹X^⊤Xw+ (X^⊤X)⁻¹X^⊤ε−w

= (X^⊤X)⁻¹X^⊤ε.

よって

( ˆw−E[ ˆw])( ˆw−E[ ˆw])^⊤ =(

(X^⊤X)⁻¹X^⊤ε) (

(X^⊤X)⁻¹X^⊤ε)_⊤

= (X^⊤X)⁻¹X^⊤εε^⊤X(X^⊤X)⁻¹. ここで，ε はノイズ項であり，

E[εε^⊤] =σ²I であることに留意すると，

E[( ˆw−E[ ˆw])( ˆw−E[ ˆw])^⊤] = (X^⊤X)⁻¹X^⊤E[εε^⊤]X(X^⊤X)⁻¹

=σ²(X^⊤X)⁻¹(X^⊤X)(X^⊤X)⁻¹

=σ²(X^⊤X)⁻¹.

演習問題 2 の解答

以下にRコードを示す：

> X <- matrix(c(1,2,2,3,4,3,2,5,3,5), nrow=5, ncol=2);

> y <- c(2,6,5,6,9)

> X

[,1] [,2]

[1,] 1 3

[2,] 2 2

[3,] 2 5

[4,] 3 3

[5,] 4 5

> y

[1] 2 6 5 6 9

> n <- nrow(X);

> d <- ncol(X);

> n [1] 5

> d [1] 2

> normX1 <- (X[,1] - mean(X[,1]))/sd(X[,1]);

> normX2 <- (X[,2] - mean(X[,2]))/sd(X[,2]);

> normX <- cbind(normX1, normX2);

> normX

normX1 normX2 [1,] -1.2278812 -0.4472136 [2,] -0.3508232 -1.1925696 [3,] -0.3508232 1.0434984 [4,] 0.5262348 -0.4472136 [5,] 1.4032928 1.0434984

> normY <- (y - mean(y))/sd(y);

> normY

[1] -1.4342743 0.1593638 -0.2390457 0.1593638 1.3545924

> XtX <- t(normX)%*%normX;

> XtY <- t(normX)%*%normY;

> invXtX <- solve(XtX);

> hat_w <- invXtX%*%XtY;

> XtX

normX1 normX2 normX1 4.000000 1.830417 normX2 1.830417 4.000000

> XtY

[,1]

normX1 3.773825 normX2 1.544176

> invXtX

normX1 normX2 normX1 0.3162162 -0.1447019 normX2 -0.1447019 0.3162162

> hat_w

[,1]

normX1 0.96989957 normX2 -0.05778621

> hat_w_orig1 <- hat_w[1]*sd(y)/sd(X[,1])

> hat_w_orig2 <- hat_w[2]*sd(y)/sd(X[,2])

> hat_w_orig0 <- mean(y) - ( hat_w_orig1*mean(X[,1])+hat_w_orig2*mean(X[,2]) )

> hat_w_orig0 [1] 0.8648649

> hat_w_orig1 [1] 2.135135

> hat_w_orig2 [1] -0.1081081

> XCovMat <- XtX/(n-1);

> XCovMat

normX1 normX2 normX1 1.0000000 0.4576043 normX2 0.4576043 1.0000000

> hat_sigma2 <- as.numeric(t(normY - normX%*%hat_w)%*%(normY - normX%*%hat_w)/(n-3));

> CovW <- invXtX*hat_sigma2;

> hat_sigma2 [1] 0.2145002

> CovW

normX1 normX2 normX1 0.06782845 -0.03103859 normX2 -0.03103859 0.06782845

> xnew <- c(1, 1);

> xnewNorm <- numeric(2);

> xnewNorm[1] <- (xnew[1] - mean(X[,1]))/sd(X[,1])

> xnewNorm[2] <- (xnew[2] - mean(X[,2]))/sd(X[,2])

> xnewNorm

[1] -1.227881 -1.937926

> EyNewNorm <- t(hat_w)%*%xnewNorm

> VyNewNorm <- t(xnewNorm)%*%CovW%*%xnewNorm

> EyNewNorm [,1]

[1,] -1.078936

> VyNewNorm [,1]

[1,] 0.2092826

> EyNew <- sd(y)*EyNewNorm + mean(y);

> VyNew <- sd(y)*sd(y)*VyNewNorm;

> EyNew

[,1]

[1,] 2.891892

> VyNew

[,1]

[1,] 1.318481

演習問題 3 の解答

S_allの計算

データが正規化されておりy¯=0であるので

S_all = (y−y)¯ ^⊤(y−y) =¯ y^⊤y.

Sregの計算

同様に，データが正規化されておりy¯=0であるので S_reg = ( ˆy−y)¯ ^⊤( ˆy−y) = ˆ¯ y^⊤y.ˆ ここで，

ˆ

y=Xwˆ =X(X^⊤X)⁻¹X^⊤y より，

S_reg = ˆy^⊤yˆ

= (X(X^⊤X)⁻¹X^⊤y)^⊤(X(X^⊤X)⁻¹X^⊤y)

=y^⊤X^⊤(X^⊤X)⁻¹X^⊤X(X^⊤X)⁻¹X^⊤y

=y^⊤X^⊤(X^⊤X)⁻¹X^⊤y.

S_resの計算

S_res= (y−y)ˆ ^⊤(y−y)ˆ

=y^⊤(I −X(X^⊤X)⁻¹X^⊤)^⊤(I−X(X^⊤X)⁻¹X^⊤)y

=y^⊤(I −X(X^⊤X)⁻¹X^⊤)y

=y^⊤y−y^⊤X(X^⊤X)⁻¹X^⊤y ただし，

y−yˆ=y−X(X^⊤X)⁻¹X^⊤y

= (I−X(X^⊤X)⁻¹X^⊤)y, ならびに

(I−X(X^⊤X)⁻¹X^⊤)^⊤(I −X(X^⊤X)⁻¹X^⊤)

=I−2X(X^⊤X)⁻¹X^⊤+X(X^⊤X)⁻¹X^⊤X(X^⊤X)⁻¹X^⊤

=I−X(X^⊤X)⁻¹X^⊤. を利用している．

以上より，

S_all =S_reg +S_res.