「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。

MENU

統計的機械学習の数理100問(02/X)

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

を用いることにする。

1. 線形回帰

 Yを目的変数、Xを説明変数として


\begin{aligned}
Y=\beta_0+\beta_1 X,\ \beta_0,\beta_1\in\mathbb{R}
\end{aligned}

という単回帰モデルを考えよう。(単)回帰モデルの場合、最小二乗法での推定が定番である。

1.1 最小二乗法によるパラメータ推定

 問題設定をしておこう。観測値として(x_1,y_1),\cdots,(x_N,y_N)が得られているとする。このとき、残差e_i


\begin{aligned}
e_i=y_{i}-\left(\beta_0+\beta_1 x_{i} \right)
\end{aligned}

で定義する。これは実際の観測値y_{i}とモデルによる推計値\beta_0+\beta_1 x_1との隔たりを表している。
 それらの総和として残差平方和を


\begin{aligned}
L=\sum_{i=1}^{N}\left\{y_{i}-\left(\beta_0+\beta_1 x_{i} \right)\right\}^2
\end{aligned}

で与える。二乗しているのは、残差の単なる総和では常に0になるからで、平方絶対和では数学的に取り扱いにくいからでもある。
 さてこの残差平方和Lを最小化するような(\beta_0,\beta_1)=(\hat{\beta}_0,\hat{\beta}_1)を求めるのが最小二乗法である。一般に2次関数が下に凸であることを踏まえれば、(\hat{\beta}_0,\hat{\beta}_1)L\beta_0および\beta_1でそれぞれ偏微分したものの連立方程式


{\displaystyle 
\begin{eqnarray}
 \left\{
    \begin{array}{l}
       \displaystyle{\frac{\partial L}{\partial \beta_0}}&=-2\displaystyle{\sum_{i=1}^{N}}\left\{y_i-\left(\beta_0+\beta_1 x_i\right)\right\}=0\\
       \displaystyle{\frac{\partial L}{\partial \beta_1}}&=-2\displaystyle{\sum_{i=1}^{N}}x_{i}\left\{y_i-\left(\beta_0+\beta_1 x_i\right)\right\}=0
    \end{array}
  \right.
\end{eqnarray}
}

の解である。したがって


{\displaystyle 
\begin{eqnarray}
 \left\{
    \begin{array}{l}
       \displaystyle{\sum_{i=1}^{N}}y_i-N\hat{\beta}_0-\hat{\beta}_1\displaystyle{\sum_{i=1}^{N}}x_i&=0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &\cdots(1)\\
       \displaystyle{\sum_{i=1}^{N}}x_{i}y_i-\hat{\beta}_0\displaystyle{\sum_{i=1}^{N}}x_{i}-\hat{\beta}_1\displaystyle{\sum_{i=1}^{N}}{x_{i}}^2&=0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &\cdots(2)
    \end{array}
  \right.
\end{eqnarray}
}

である。
 まず


\begin{aligned}
\bar{x}&=\displaystyle{\frac{1}{N}\sum_{i=1}^{N}x_i},\\
\bar{y}&=\displaystyle{\frac{1}{N}\sum_{i=1}^{N}y_i}
\end{aligned}

とおくと、(1)より


\begin{aligned}
&N\hat{\beta}_0=\displaystyle{\sum_{i=1}^{N}}y_i-\hat{\beta}_1\displaystyle{\sum_{i=1}^{N}}x_i\\
\Leftrightarrow&\hat{\beta}_0=\bar{y}-\hat{\beta}_1\bar{x}
\end{aligned}

である。
 次に(1)\times\left(\displaystyle{\sum_{i=1}^{N}} x_i\right)-(2)\times Nより


\begin{aligned}
\left(\displaystyle{\sum_{i=1}^{N}} x_i\right)\left(\displaystyle{\sum_{i=1}^{N}} y_i\right)-N\left(\displaystyle{\sum_{i=1}^{N}} x_i y_i\right)+\hat{\beta}_1\left\{N\displaystyle{\sum_{i=1}^{N}} {x_i}^2 -\left(\displaystyle{\sum_{i=1}^{N}} x_i \right)^2\right\}=0\\
\therefore\ \ \hat{\beta}_1=\displaystyle{\frac{N\left(\displaystyle{\sum_{i=1}^{N}} x_i y_i \right)-\left(\displaystyle{\sum_{i=1}^{N}} x_i\right)\left(\displaystyle{\sum_{i=1}^{N}} y_i\right)}{N\left(\displaystyle{\sum_{i=1}^{N}} {x_i}^2\right)-\left(\displaystyle{\sum_{i=1}^{N}} x_i \right)^2}}
\end{aligned}

分母分子をともにNで割ることにより、


\begin{aligned}
\hat{\beta}_1&=\displaystyle{\frac{\left(\displaystyle{\sum_{i=1}^{N}} x_i y_i \right)-\bar{x}\left(\displaystyle{\sum_{i=1}^{N}} y_i\right)}{\left(\displaystyle{\sum_{i=1}^{N}} {x_i}^2\right)-\bar{x}\displaystyle{\sum_{i=1}^{N}} x_i}}\\
&=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{N}}\left(x_i y_i-\bar{x}y_i\right)}{\displaystyle{\sum_{i=1}^{N}}\left({x_i}^2-\bar{x}x_i\right)}}\\
&=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{N}}y_i\left(x_i-\bar{x}\right)}{\displaystyle{\sum_{i=1}^{N}}x_i\left(x_i-\bar{x}\right)}}
\end{aligned}

ここで


\begin{aligned}
\bar{x}\displaystyle{\sum_{i=1}^{N}\left(x_i-\bar{x}\right)}=0,\ \ \bar{y}\displaystyle{\sum_{i=1}^{N}\left(x_i-\bar{x}\right)}=0
\end{aligned}

であるから、前者を分母に後者を分子にそれぞれ加算することで


\begin{aligned}
\hat{\beta}_1=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{N}}\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}{\displaystyle{\sum_{i=1}^{N}}\left(x_i-\bar{x}\right)^2}}
\end{aligned}

が成り立つ。以上をまとめれば、


\begin{aligned}
\hat{\beta}_1&=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{N}}\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}{\displaystyle{\sum_{i=1}^{N}}\left(x_i-\bar{x}\right)^2}},\\
\hat{\beta}_0&=\bar{y}-\hat{\beta}_1\bar{x}
\end{aligned}

である。

1.2 中心化

 前節で求めた(\hat{\beta}_0,\hat{\beta}_1)について、


\begin{aligned}
x_{i}^{\prime}=x_{i}-\bar{x},\ \ y_{i}^{\prime}=y_{i}-\bar{y}
\end{aligned}

と置き換えるとx_{i}=x_{i}^{\prime}+\bar{x},\ \ y_{i}=y_{i}^{\prime}+\bar{y}であり


\begin{aligned}
\bar{x}_{i}^{\prime}=\bar{y}_{i}^{\prime}=0
\end{aligned}

であるから


\begin{aligned}
\hat{\beta}_1^{\prime}&=\displaystyle{\frac{\displaystyle{\sum_{i=1}^{N}}x_{i}^{\prime}y_{i}^{\prime}}{\displaystyle{\sum_{i=1}^{N}}{x_{i}^{\prime}}^2}},\\
\hat{\beta}_0^{\prime}&=\bar{y}_{i}^{\prime}-\hat{\beta}_1^{\prime}\bar{x}_{i}^{\prime}=0
\end{aligned}

が成り立つ。
 これは座標平面上に(x_i,y_i)をプロットしたものを考えたときに、それらをx方向に\bar{x}y方向に\bar{y}にだけ平行移動したものである。すなわち直線


\begin{aligned}
y&=\hat{\beta}_0+\hat{\beta}_1x,\\
y&=\hat{\beta}_0^{\prime}+\hat{\beta}_1^{\prime}x=\hat{\beta}_1^{\prime}x
\end{aligned}

について、後者は前者の傾きをそのままに原点を通るように平行移動したものを表すことになる。

1.3 RおよびPythonによる実装

 入力データは以下のCSVファイルを用いた:

https://drive.google.com/file/d/10iv9jZYqvYpFJx2wiRPxJi2Lyiu-essp/view?usp=sharing

1.3.1 Rの実装
##########################################################
### 最小二乗法により単回帰モデルのパラメータを推定する ###
##########################################################

### 真の値:
### beta_0 = 4.2
### beta_1 = 0.5


# 初期設定
options(stringsAsFactors = F)


# ディレクトリの設定
dir.input <- "...適当な入力ファイルのあるディレクトリを指定..."
dir.output <- "...適当な出力先を指定…"


# パラメータ推定の関数
EstParamLM <- function(vc.x,vc.y){
  x.mean <- mean(vc.x)
  y.mean <- mean(vc.y)
  
  # β_1の推定
  beta.1 <- (vc.x-x.mean) %*% (vc.y-y.mean)
  beta.1 <- beta.1/((vc.x-x.mean) %*% (vc.x-x.mean))
  
  # β_0の推定
  beta.0 <- y.mean - beta.1 * x.mean
  
  return(c(beta.0, beta.1))
}


# 実際に推定する
df.input <- read.csv(paste0(dir.input,"01_01_TestData.csv"),header = T)

vc.start <- Sys.time()
vc.EstParam <- EstParamLM(df.input[,"x"],df.input[,"y"])
vc.end <- Sys.time() - vc.start

# 推定結果と計算時間を出力
vc.AbsError <- vc.EstParam-c(4.2,0.5) # 絶対誤差
vc.RelError <- vc.AbsError/c(4.2,0.5) # 絶対誤差

vc.out <- c("R",vc.EstParam,vc.AbsError,vc.RelError,as.numeric(vc.end))
df.output <- as.data.frame(t(vc.out))

colnames(df.output) <- c("言語","EstBeta0","EstBeta1","AbsErrorBeta0","AbsErrorBeta1","RelErrorBeta0","RelErrorBeta1","CalcTime[s]")

write.csv(x = df.output, file = paste0(dir.output,"01_01_Output_R.csv"),row.names = F)
1.3.2 Pythonの実装
import csv
import numpy as np
import time

def EstParamLM(x,y):
    x_bar, y_bar = np.mean(x), np.mean(y)
    beta1 = np.dot(x - x_bar, y - y_bar)/np.linalg.norm(x - x_bar) ** 2
    beta0 = y_bar - beta1 * x_bar
    return [beta0, beta1]

# 読み込みファイルの整理
dir_input <- '...適当な入力ファイルのあるディレクトリを指定...'
dir_output <- '...適当な出力先を指定…'

file_input = open(dir_input + '01_01_TestData.csv', 'r', encoding='utf-8', errors='', newline='')

#リスト形式
f = csv.reader(file_input, delimiter=',', doublequote=True, lineterminator='\r\n', quotechar='"', skipinitialspace=True)
header = next(f)  # ヘッダーを読み飛ばす

lst_z = np.asarray([row for row in f], dtype=np.float32)

lst_x = np.asarray([row[0] for row in lst_z], dtype=np.float32)
lst_y = np.asarray([row[1] for row in lst_z], dtype=np.float32)

start = time.time()
beta0,beta1 = EstParamLM(lst_x,lst_y)
elapsed_time = time.time() - start

print(beta0),print(beta1)

with open(dir_output + '01_01_Output_Python.csv', 'w') as f:
    writer = csv.writer(f)
    writer.writerow(['言語','EstBeta0','EstBeta1','AbsErrorBeta0','AbsErrorBeta1','RelErrorBeta0','RelErrorBeta1','CalcTime[s]'])
    writer.writerow(['Python',beta0, beta1, beta0-4.2, beta1-0.5,(beta0-4.2)/4.2, (beta1-0.5)/0.5,format(elapsed_time)])

1.3.3 計測結果

 上記をそれぞれ実行したときの結果を以下に掲げる:

言語
\hat{\beta}_0
\hat{\beta}_1
絶対誤差(\beta_0)
絶対誤差(\beta_1)
相対誤差(\beta_0)
相対誤差(\beta_1)
計測時間
R 4.20000000006256 0.499999999921655 6.26\times10^{-11} -7.83\times10^{-11} 1.49\times10^{-11} -1.57\times10^{-10} 0.001185894
Python 4.200000014 0.499999947 1.36\times 10^{-8} -5.32\times 10^{-8} 3.23\times 10^{-9} -1.06\times 10^{-7} 0
1.3.3.1 PCスペック情報
エディション Windows 10 Home
バージョン 20H2
プロセッサ Intel(R) Core(TM) i5-1035G4 CPU @ 1.10GHz 1.50 GHz
実装 RAM 8.00 GB
システムの種類 64 ビット オペレーティング システム、x64 ベース プロセッサ
R バージョン 3.6.3 (2020-02-29)
RStudio バージョン 1.2.5033
Python バージョン 3.7.9 64-bit
Spyder バージョン 5.1.5 None
プライバシーポリシー お問い合わせ