この資料では、線形回帰モデルの回帰係数をPython、R、Juliaの3つの言語で推定する方法について紹介する。どの言語でもライブラリを使う方法と、NumPyやそれに相当するような機能のみを用いたやり方を紹介している。
筆者はJuliaのDataFrames.jlや行列計算に関して調べ物が必要だったので、以下のサイトに大変お世話になった。ありがたい。ほかの2つの言語については、別の機会に参考になりそうなものをまとめる。
また、線形回帰に関する話題は線形回帰モデルに関する説明をしているサイト、の欄をみてほしい。説明をみごとにまるごとブン投げてしまっているのでリッジ回帰あたりからはちゃんと書き直したいよねの気持ち。
ここでは、線形回帰モデルの説明をするとみせかけて参考となるようなサイトを紹介するにとどめる。線形回帰モデルについての説明は2022年4月時点でも素晴らしいものが星の数ほどある状態で、いい参考文献を紹介する方が筆者の曖昧な理解から繰り出される怪文書よりもよほどいいと思ったからそうする。
Rによる回帰係数の推定Rによる回帰係数の推定は、よくlm関数によって行われている。そこで、本資料では下記の2種類の方法で回帰係数を推定してその結果を比較する。
lm関数による推定Rの行列計算による推定lm関数による推定# irisを直接用いる
model <- lm("Sepal.Width ~ Sepal.Length + Petal.Length", data = iris)
model$coef## (Intercept) Sepal.Length Petal.Length
## 1.0380691 0.5611860 -0.3352667Rの行列計算による推定my_coefは、model$coefと誤差などの情報以外は同じ結果であることがわかる。
# データの処理はなるべくtidyverseを用いることにした。
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.1 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()X0 <- iris %>% select(Sepal.Length, Petal.Length)
intercept <- rep(1, dim(X0)[1])
X <- cbind(intercept, X0) %>% as.matrix()
y <- iris %>% select(Sepal.Width) %>% as.matrix()
my_coef <- solve(t(X)%*%X)%*%t(X)%*%y
my_coef## Sepal.Width
## intercept 1.0380691
## Sepal.Length 0.5611860
## Petal.Length -0.3352667Pythonによる推定Pythonによる回帰係数の推定は、ライブラリのsklearnかstatsmodelsによって行われている。この2つのライブラリによる推定結果を比較すると、statsmodelsのほうが推定に関する情報を多く記録してくれるという印象がある。そのため、本資料では下記の2種類の方法で回帰係数を推定してその結果を比較する。
statsmodelsによる推定NumPyの行列計算による推定statsmodelsによる推定推定結果model.paramsは、my_coefやmodel$coefと誤差を除けば同じになっている。
# インポート
import statsmodels.api as sm
# データirisをインポートする
iris = sm.datasets.get_rdataset("iris", "datasets").data
# 説明変数の行列Xと被説明変数のベクトルyを作る
X = iris[["Sepal.Length", "Petal.Length"]]
y = iris["Sepal.Width"]
# sm.add_constant(X)は、定数項部分の1をXに加えたもの。
model = sm.OLS(y, sm.add_constant(X)).fit()
# 係数をとってきたい場合、paramsメソッドを用いる。
model.params## const 1.038069
## Sepal.Length 0.561186
## Petal.Length -0.335267
## dtype: float64NumPyの行列計算による推定推定結果coefは、いままで行ってきた推定結果に一致しているように見える。
# インポート
import numpy as np
import statsmodels.api as sm
# データirisをインポートする
iris = sm.datasets.get_rdataset("iris", "datasets").data
# 説明変数の行列Xと被説明変数のベクトルyを作る
# sm.add_constant(X)と同じことをするためにiris["const"]を追加
iris["const"] = 1
X = iris[["const", "Sepal.Length", "Petal.Length"]].values
y = iris["Sepal.Width"].values
# 演算子@で内積をとり、np.linalg.solveで引数1の逆行列を引数2にかけたものを計算できる。
# 下の書き方は速度面に配慮した記法で、次の記事を参考にした。
# https://qiita.com/fujiisoup/items/e7f703fc57e2dfc441ad
coef = np.linalg.solve(X.T @ X, X.T @ y)
coef## array([ 1.03806906, 0.56118597, -0.33526674])Juliaによる推定Juliaによる回帰係数の推定は、ライブラリのGLMが用意しているlm関数によって行われることが多いようだ。そこで、本資料では下記の2種類の方法で回帰係数を推定してその結果を比較する。
GLMのlm関数による推定LinearAlgebraの行列計算による推定GLMによる推定値はいままで推定してきた結果とほとんど変わらないことがわかる。
# インポート
using RDatasets
using GLM
# RDataSetsからirisをロードする
df = dataset("datasets", "iris");
model = lm(@formula(SepalWidth ~ SepalLength + PetalLength), df);
GLM.coef(model)## 3-element Vector{Float64}:
## 1.0380690570704476
## 0.5611859710980439
## -0.33526674157885833LinearAlgebraを使った行列計算による推定LinearAlgebraによる回帰係数の推定を行い、ほかの値と比較して差がないことを確認している。
# インポート
using RDatasets
using DataFrames
using LinearAlgebra
# RDataSetsからirisをロードする
df = dataset("datasets", "iris");
df[!, :const] = [1 for i = 1:size(df)[1]];
X = Matrix(df[:, [:const, :SepalLength, :PetalLength]]);
y = df[:, :SepalWidth];
# 行列計算を行う。
coef = inv(transpose(X) * X)* transpose(X) * y;
coef## 3-element Vector{Float64}:
## 1.038069057070371
## 0.5611859710980548
## -0.33526674157886044