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

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

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

作者:讓, 鈴木
共立出版

Amazon

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

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回

power-of-awareness.com

1. 線形回帰

1.4 重回帰

　説明変数が $p\gt1$ 個ある場合の回帰を重回帰と呼ぶがそれを検討する。

$\begin{aligned} \boldsymbol{y}=\begin{bmatrix} y_{1} \\ \vdots \\ y_{p} \end{bmatrix}, X=\begin{bmatrix} 1 & x_{11} & \cdots & x_{1N} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{p1} & \cdots & x_{pN} \end{bmatrix}, \boldsymbol{\beta}=\begin{bmatrix} \beta_0 \\ \vdots \\ \beta_p \end{bmatrix} \end{aligned}$

とすれば、

$\begin{aligned} L=\|\boldsymbol{y}-X\boldsymbol{\beta}\|^2 \end{aligned}$

と書け、

$\begin{aligned} \nabla{L}=\begin{bmatrix} \displaystyle{\frac{\partial }{\partial \beta_0}} \\ \vdots \\ \displaystyle{\frac{\partial }{\partial \beta_p}} \end{bmatrix}=-2\ ^tX\left(\boldsymbol{y} -X\boldsymbol{\beta}\right) \end{aligned}$

が成立する。
　単回帰のときと同様に、 $^{t}XX$ が逆行列をもつならば

$\begin{aligned} \hat{\boldsymbol{\beta}}=\left(^{t}XX \right)^{-1}\ ^{t}X\boldsymbol{y} \end{aligned}$

が成り立つ。

1.4.2 RおよびPythonによる実装

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

https://drive.google.com/file/d/1VqaXF2xcBme87bLrGXsQUHgD5t1QUVFa/view?usp=sharing

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

# モデル
# y_i=beta_0+x_1*beta_1+x_2*beta_2+x_3*beta_3+eps_i

### 真の値：
# beta_0 = -0.5
# beta_1 = 3.1
# beta_2 = -2.1
# beta_3 = 1.1
# eps~N(0,1)


# 初期設定
options(stringsAsFactors = F)


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


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

vc.ObjVar <- df.input[,ncol(df.input)]
mt.ExpVars <- as.matrix(df.input[,-ncol(df.input)])

vc.start <- Sys.time()
vc.beta <- solve(t(mt.ExpVars)%*%mt.ExpVars)%*%t(mt.ExpVars)%*%vc.ObjVar
vc.end <- Sys.time() - vc.start

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

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

colnames(df.output) <- c("言語","EstBeta0","EstBeta1","EstBeta2","EstBeta3","AbsErrorBeta0","AbsErrorBeta1","AbsErrorBeta2","AbsErrorBeta3","RelErrorBeta0","RelErrorBeta1","RelErrorBeta2","RelErrorBeta3","CalcTime[s]")

write.csv(x = df.output, file = paste0(dir.output,"01_02_Output_R.csv"),row.names = F)

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

# モデル
# y_i=beta_0+x_1*beta_1+x_2*beta_2+x_3*beta_3+eps_i

### 真の値：
# beta_0 = -0.5
# beta_1 = 3.1
# beta_2 = -2.1
# beta_3 = 1.1
# eps~N(0,1)

import csv
import itertools
import numpy as np
import pandas as pd
import time

# ディレクトリの設定
dir_input = "...入力ファイルのある適当なパス..."
dir_output = "...適当な出力先のパス..."

# データを読み込み
file_input = open(dir_input + '01_02_TestData.csv', 'r', encoding='utf-8', errors='', newline='')

#リスト形式
f = pd.read_csv(file_input)

# xとyに分離する
f_x = f.drop('y',axis = 1)
f_y = f['y']

# 推定する
start = time.time()
betas = np.linalg.inv(f_x.T@f_x)@f_x.T@f_y
elapsed_time = time.time() - start

# 出力
output = list(itertools.chain.from_iterable([[beta for beta in betas], [beta for beta in betas-[-0.5,3.1,-2.1,1.1]],[beta for beta in (betas-[-0.5,3.1,-2.1,1.1])/[-0.5,3.1,-2.1,1.1]]]))
output.insert(0, 'Python')
output.insert(len(output),elapsed_time)

with open(dir_output + '01_02_Output_Python.csv', 'w', newline="") as f:
    writer = csv.writer(f)
    writer.writerow(['言語','EstBeta0','EstBeta1','EstBeta2','EstBeta3','AbsErrorBeta0','AbsErrorBeta1','AbsErrorBeta2','AbsErrorBeta3','RelErrorBeta0','RelErrorBeta1','RelErrorBeta2','RelErrorBeta3','CalcTime[s]'])
    writer.writerow(output)

1.4.2.1　スペック情報

エディション	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