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

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年12月10日改訂)。

MENU

統計的機械学習の数理100問(03/20)

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

を用いることにする。

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

1.4.2.2 計算結果

言語
\hat{\beta}_0
\hat{\beta}_1
\hat{\beta}_2
\hat{\beta}_3
絶対誤差(\beta_0)
絶対誤差(\beta_1)
絶対誤差(\beta_2)
絶対誤差(\beta_3)
相対誤差(\beta_0)
相対誤差(\beta_1)
相対誤差(\beta_2)
相対誤差(\beta_3)
計測時間
R -0.499999999528499 3.09999999979884 -2.10000000056355 1.09999999951428 4.72\times 10^{-10} 2.01\times 10^{-10} 5.64\times 10^{-10} 4.86\times 10^{-10} -9.43\times 10^{-10} -6.49\times 10^{-11} 2.68\times 10^{-10} -4.42\times10^{-10} 0.001991987
Python -0.499999999528501 3.09999999979884 -2.10000000056354 1.09999999951427 4.71\times10^{-10} 2.01\times10^{-10} 5.64\times10^{-10} 4.86\times10^{-10} -9.43\times10^{-10} -6.49\times10^{-11} 2.68\times10^{-10} -4.42\times10^{-10} 0.002031326
プライバシーポリシー お問い合わせ