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

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

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

作者:讓, 鈴木
共立出版

Amazon

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

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回

power-of-awareness.com

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

次回