Screaming Loud

日々是精進

scala(breeze)で重回帰を実装

重回帰とは

重回帰は色んな所に転がっているので今更書かなくてもって感じですが、簡単に。


ある地域の住宅価格を求めたいとなった時に、ランダムで価格を当てるよりも何かの影響を受けて価格が変動すると考えるのが自然です。

 {住宅価格 = \beta_0 + \beta_1 *  犯罪率 + \beta_2 * 宅地の割合 + ...  \beta_n * 低所得者の割合 - (1)}

そこで上のような住宅価格は色んな要素に依存していると考え、各要素がどれくらい住宅価格を求めるのに影響があるのかを求めるのが重回帰です。

住宅価格に犯罪率の寄与が少なければ、{\beta_1}の値は小さくなるでしょうし、大きければその逆になります。

回帰を計算すると、各要素の係数であるβらが求まります。
βらが求まると、上記 (1) の式が埋まるので、新しいデータを各要素に当てはめると、そのデータに対する住宅価格が推定できるようになります。

この住宅価格を目的変数、犯罪率などの要素を説明変数と呼びます。

詳しくは、R - 重回帰モデルの理論と実装 -なぜ正則化が必要か- - Qiitaなどで式の導出は詳しく書かれています。
気になる方は、ぜひ導出も追ってみてください。

実装

Bostonの住宅の価格のCSVです。
Sense · Enterprise Data Science Platform

Python

上述のページで書かれているPythonコードの引用が以下になります。

import pandas as pd
import numpy as np
from numpy import linalg as la

df = pd.read_csv("Boston.csv", index_col=0)
y = df.iloc[:,  13].values
df = (df - df.mean())/df.std() # 基準化
X = df.iloc[:, :13].values
X = np.column_stack((np.ones(len(X)),X)) 
beta = np.dot(np.dot(la.inv(np.dot(X.T, X)),X.T), y)
print beta
R

Rのlm関数で書いたものの引用です

library(MASS)
X <- as.matrix(Boston[, -14])
y <- as.vector(Boston[,  14])
X <- scale(X) #基準化
model <- lm(y ~ X)
print(model)
scala

この式をScalaで書き直します。

def regression(file: File): DenseVector[Double] = {
  val mat = csvread(file, skipLines = 1)
  val y = mat(::, mat.cols - 1) // 目的変数
  val X = DenseMatrix.horzcat(
    DenseMatrix.tabulate(mat.rows, 1) { case _ => 1.0 }, // 要素がない列
    mat(::, 1 to mat.cols - 2) // 要素がある列
  ) // 説明変数の行列
  (inv(X.t * X) * X.t) * y
}

出力

DenseVector(36.45948838509085, -0.10801135783679613, 0.04642045836688007, 0.020558626367099816, 2.6867338193447603, -17.76661122830122, 3.809865206809127, 6.92224640347558E-4, -1.4755668456002482, 0.3060494789851876, -0.012334593916574491, -0.9527472317074077, 0.009311683273793751, -0.5247583778554944)
解説

上の式は

 {住宅価格 = \beta_0 + \beta_1 *  犯罪率 + \beta_2 * 宅地の割合 + ...  \beta_n * 低所得者の割合 - (1)}

Index,犯罪率,宅地の割合,...,低所得者の割合,住宅価格
1,0.5, 0.8, ... , 0.2, 5003,
2,0.53, 0.82, ... , 0.12, 600
3,...

のようなCSVを受け取ります。

まず、1行目で、CSVをロードして行列クラスにぶち込みます。

yは目的変数です。
上のcsvをそのまま行列に突っ込むと、最終列が住宅価格になるので、最終行を目的変数のベクトルとします。

Xは各要素群ですが、そのままだと{\beta_0}の列が作られません。
なので、tabulateで作成して、横で合成しています。
ここは見やすく、fillでVectorを作ってtoDenseMatrixでもよいです。

あとは、最小二乗法の解を求めるだけです。
{\beta = (X^TX)^{-1}X^Ty}
細かい導出式は、今回は求めませんので、他のサイトを見てみるとよいと思います。


今回は基準化を省いていますが、基準化を行うと各要素が同じスケールになるので、上の出力に出た係数を比較するだけでどの要素が目的変数に影響を及ぼすかがわかります。

ただその場合、係数のスケールが変わってしまっているので、目的変数の導出はスケールを戻さないと求められなくなります。


以下に動く実装を置いています。
regression/Regression.scala at master · moc-yuto/regression · GitHub