統計的機械学習の数理100問（07/20） - 「大人の教養・知識・気付き」を伸ばすブログ

　いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回
1. 線形回帰
- 1.8.　信頼区間と予測区間
  - 1.8.1　シミュレーション(R)
  - 1.8.2　シミュレーション(Python)
次回

前回

power-of-awareness.com

1. 線形回帰

1.8.　信頼区間と予測区間

　 $\boldsymbol{\beta}\in\mathbb{R}^{p+1}$ の推定量 $\hat{\boldsymbol{\beta}}$ の信頼区間を考える。 $SE(\cdot)$ を標準誤差を表す演算子、 $n$ を標本サイズだとして

$\begin{aligned} t_i=\displaystyle{\frac{\hat{\beta}_i-\beta_i}{SE(\hat{\beta}_i)}}\sim t(n-p-1) \end{aligned}$

が成り立つから、 $\beta_i$ の確率密度関数を $f(\cdot)$ として

$\begin{aligned} \displaystyle{\frac{\alpha}{2}}=\displaystyle{\int_{t}^{\infty}f(u)}du \end{aligned}$

を満たすような $t$ を $t_{n-p-1}\left(\displaystyle{\frac{\alpha}{2}}\right)$ とおけば $\hat{\boldsymbol{\beta}}$ の信頼区間として

$\begin{aligned} \beta_i\in\left[\hat{\beta}_i-t_{n-p-1}\left(\displaystyle{\frac{\alpha}{2}}\right) SE(\hat{\beta}_i),\hat{\beta}_i+t_{n-p-1}\left(\displaystyle{\frac{\alpha}{2}}\right) SE(\hat{\beta}_i)\right] \end{aligned}$

が得られる。
　推定に用いた標本 $\boldsymbol{X}_1,\cdots,\boldsymbol{X}_n$ とは別に $\boldsymbol{X}_{*}\in\mathbb{R}^{p+1}$ を用いて目的変数を予測するのに $\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}$ を用いたとする。この予測の信頼区間を考える。このとき

$\begin{aligned} E[\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}]=\boldsymbol{X}_{*} E[\hat{\boldsymbol{\beta}}] \end{aligned}$

が成り立つから、その分散は、 $V[\varepsilon_i]=\sigma^2,i=1,2,\cdots,n$ として

$\begin{aligned} V[\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}]&=E[{}^{t}(\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-E[\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}])(\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-E[\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}])]\\ &=E[{}^{t}(\hat{\boldsymbol{\beta}}-E[\hat{\boldsymbol{\beta}}]){}^{t}\boldsymbol{X}_{*}\boldsymbol{X}_{*}(\hat{\boldsymbol{\beta}}-E[\hat{\boldsymbol{\beta}}])]\\ &=E[{}^{t}(\hat{\boldsymbol{\beta}}-E[\hat{\boldsymbol{\beta}}]){}^{t}\boldsymbol{X}_{*}\boldsymbol{X}_{*}(\hat{\boldsymbol{\beta}}-E[\hat{\boldsymbol{\beta}}])]\\ &=\boldsymbol{X}_{*}\mathbb{V}[\hat{\boldsymbol{\beta}}]{}^{t}\boldsymbol{X}_{*}\\ &=\sigma^2 \boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X} \boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*} \end{aligned}$

を得る。

$\begin{aligned} \hat{\sigma}&=\sqrt{\displaystyle{\frac{RSS}{n-p-1}}},\\ SE(\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}})&=\hat{\sigma}\sqrt{\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X} \boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}} \end{aligned}$

とおくとき、

$\begin{aligned} C&=\displaystyle{\frac{\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-\boldsymbol{X}_{*}\boldsymbol{\beta}}{SE(\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}})}}\\ &=\displaystyle{\frac{\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-\boldsymbol{X}_{*}\boldsymbol{\beta}}{\hat{\sigma}\sqrt{\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X} \boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}}}}\\ &=\displaystyle{\frac{\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-\boldsymbol{X}_{*}\boldsymbol{\beta}}{\sigma\sqrt{\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X} \boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}}}\left/\sqrt{\displaystyle{\frac{RSS}{\sigma^2}/(n-p-1)}}\right.} \end{aligned}$

が自由度 $n-p-1$ の $t$ 分布に従うことが知られている。
　他方で $\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}$ で推定した値に、実際には誤差項 $\boldsymbol{\varepsilon}$ が加味されるのでその値を考慮した場合に $\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}$ と $Y_{*}=\boldsymbol{X}_{*}\boldsymbol{\beta}+\boldsymbol{\varepsilon}$ との差の分布を評価する必要がある。この場合の信頼区間は予測区間と呼ばれる。

$\begin{aligned} V[\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-(\boldsymbol{X}_{*}\boldsymbol{\beta}+\boldsymbol{\varepsilon})]&=V[\boldsymbol{X}_{*}\left(\hat{\boldsymbol{\beta}}\boldsymbol{\beta}\right)]+V[\boldsymbol{\varepsilon})]\\ &=\sigma^2 \boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X} \boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}+\sigma^2 \end{aligned}$

であり、以上と同様にして

$\begin{aligned} P&=\displaystyle{\frac{\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-Y_{*}}{SE(\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}})}}\\ &=\displaystyle{\frac{\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-Y_{*}}{\hat{\sigma}\sqrt{\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X} \boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}}}}\\ &=\displaystyle{\frac{\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-Y_{*}}{\sigma\left(1+\sqrt{\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X} \boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}}\right)}\left/\sqrt{\displaystyle{\frac{RSS}{\sigma^2}\left/(n-p-1)\right.}}\right.}\sim t(n-p-1) \end{aligned}$

を得る。したがって棄却域を信頼水準 $\alpha$ を用いて

$\begin{aligned} \boldsymbol{X}_{*}\boldsymbol{\beta}\in&\left[\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-t_{n-p-1}\left(\displaystyle{\frac{\alpha}{2}}\right)\hat{\sigma}\sqrt{\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}},\right.\\&\left.\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}+t_{n-p-1}\left(\displaystyle{\frac{\alpha}{2}}\right)\hat{\sigma}\sqrt{\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}}\right],\\ \hat{Y}_{*}&\in\left[\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}-t_{n-p-1}\left(\displaystyle{\frac{\alpha}{2}}\right)\hat{\sigma}\sqrt{1+\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}}\right.,\\&\left.\boldsymbol{X}_{*}\hat{\boldsymbol{\beta}}+t_{n-p-1}\left(\displaystyle{\frac{\alpha}{2}}\right)\hat{\sigma}\sqrt{1+\boldsymbol{X}_{*}\left({}^{t}\boldsymbol{X}\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}_{*}}\right] \end{aligned}$

を得る。

1.8.1　シミュレーション(R)

############################################
### 信頼区間と予測区間のシミュレーション ###
############################################

library("ggplot2")
library("reshape2")

set.seed(1) # 再現性の確保

### 標本の作成
num_sample <- 1000
p <- 1

# 標本ベクトルの作成
mt_sample <- matrix(rnorm(num_sample * p), ncol = p)
mt_sample <- cbind(rep(1,num_sample),mt_sample)

# 係数ベクトルの作成
vc_beta <- rnorm(p + 1, 0, 1)

# 誤差項の作成
vc_epsilon <- rnorm(num_sample)

# 目的変数の作成
mt_obj <- mt_sample %*% vc_beta + vc_epsilon

### 推定値の計算
mt_inv <- solve(t(mt_sample)%*%mt_sample)
beta_hat <- mt_inv %*% t(mt_sample) %*% mt_obj

# RSSを計算
num_RSS <- sum((mt_obj-mt_sample %*% beta_hat)^2)
num_RSE <- sqrt(num_RSS/(num_sample - p - 1))
num_alpha <- 0.05 # 信頼水準:5%

conf_intvl <- function(x){
  x <- cbind(1,x)
  vc_range <- qt(df = num_sample - p - 1,1 - 0.5 * num_alpha) * num_RSE * sqrt(x %*% mt_inv %*% t(x))
  
  return(list(lower = x%*%beta_hat - vc_range, upper = x%*%beta_hat + vc_range))
}

pred_intvl <- function(x){
  x <- cbind(1,x)
  vc_range <- qt(df = num_sample - p - 1,1 - 0.5 * num_alpha) * num_RSE * sqrt(1 + x %*% mt_inv %*% t(x))
  
  return(list(lower = x%*%beta_hat - vc_range, upper = x%*%beta_hat + vc_range))
}

### 信頼区間と予測区間を図示
vc_seq <- seq(-8,8,0.1)

# 信頼区間
vc_seq_conf_lower <- NULL
vc_seq_conf_upper <- NULL

for(vc in vc_seq){
  vc_seq_conf_lower <- c(vc_seq_conf_lower, conf_intvl(vc)$lower)
  vc_seq_conf_upper <- c(vc_seq_conf_upper, conf_intvl(vc)$upper)  
}

# 予測区間
vc_seq_pred_lower <- NULL
vc_seq_pred_upper <- NULL

for(vc in vc_seq){
  vc_seq_pred_lower <- c(vc_seq_pred_lower, pred_intvl(vc)$lower)
  vc_seq_pred_upper <- c(vc_seq_pred_upper, pred_intvl(vc)$upper)  
}

# グラフ作成
vc_x_lim <- c(min(vc_seq),max(vc_seq))
vc_y_lim <- c(min(vc_seq_conf_lower,vc_seq_pred_lower),max(vc_seq_conf_upper,vc_seq_pred_upper)) 


df <- data.frame(x = vc_seq, cnf_lw = vc_seq_conf_lower, cnf_up = vc_seq_conf_upper,
                 prd_lw = vc_seq_pred_lower, prd_up = vc_seq_pred_upper
                 )
df_melt <- melt(df,id.vars = "x")

# グラフ表示
g <- ggplot(data = df_melt, aes(x = x, y= value, color = variable))
g <- g + geom_line() + geom_abline(intercept = beta_hat[1], slope = beta_hat[2])
g <- g + scale_color_manual(values = c("blue", "blue","red","red"))
g <- g + ggtitle("信頼区間と予測区間") + 
  labs(x = "x",y = "y") + theme(plot.title = element_text(hjust = 0.5),legend.position = "none",
                                legend.title=element_text(size = 7),
                                legend.text=element_text(size = 7))
g <- g + xlim(vc_x_lim) + ylim(vc_y_lim)
plot(g)

1.8.2　シミュレーション(Python)

import numpy as np
import random
from scipy import stats
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats

# 標本の組成
num_sample = 1000; p=1
X = np.random.randn(num_sample, p)
X = np.insert(X, 0, 1, axis = 1)
beta = np.array([1,1])
epsilon = np.random.randn(num_sample)
y = X@beta + epsilon

# 信頼区間・予測区間の生成
U = np.linalg.inv(X.T@X)
beta_hat = U@X.T@y
RSS = np.linalg.norm(y - X@beta_hat)**2
RSE = np.sqrt(RSS / (num_sample-p-1))
alpha=0.05

# 信頼区間
def cnf_intvl(x): 
    x = np.array([1,x])
    # stats.t.ppf(0.975,df=N-p-1) # 累積確率が1-0.5 * alphaとなる点
    range = stats.t.ppf(1-0.5*alpha, df = num_sample-p-1) * RSE * np.sqrt(x@U@x.T)
    lower = x@beta_hat - range
    upper = x@beta_hat + range
    return ([lower, upper])

# 予測区間
def prd_intvl(x): 
    x = np.array([1,x])
    # stats.t.ppf(0.975,df=N-p-1) # 累積確率が1-alpha/2となる点
    range = stats.t.ppf(1-0.5*alpha,df = num_sample-p-1) * RSE * np.sqrt(1+x@U@x.T)
    lower = x@beta_hat - range
    upper = x@beta_hat + range
    return ([lower,upper])

# 例
print(stats.t.ppf(1-0.5*alpha,df=1))  # 確率pに対応する点

x_seq = np.arange(-10,10,0.1) # x軸の表示区間

# 信頼区間
lower_seq1 = []; upper_seq1 = []
for i in range(len(x_seq)):
    lower_seq1.append(cnf_intvl(x_seq[i])[0])
    upper_seq1.append(cnf_intvl(x_seq[i])[1])
# 予測区間
lower_seq2 = []; upper_seq2 = []
for i in range(len(x_seq)):
    lower_seq2.append(prd_intvl(x_seq[i])[0])
    upper_seq2.append(prd_intvl(x_seq[i])[1])

yy = beta_hat[0] + beta_hat[1] * x_seq

plt.xlim(np.min(x_seq),np.max(x_seq))
plt.ylim(np.min(lower_seq1),np.max(upper_seq1))

plt.plot(x_seq,yy,c="black")
plt.plot(x_seq,lower_seq1,c="blue")
plt.plot(x_seq,upper_seq1,c="red")
plt.plot(x_seq,lower_seq2,c="blue",linestyle="dashed")
plt.plot(x_seq,upper_seq2,c="red",linestyle="dashed")

plt.xlabel("x")
plt.ylabel("y")