数式中の記号 | コード中の変数 | 概要 |
---|---|---|
y | 目的変数で、実数値をとる。 | |
y[i-1] | 番目の標本が持つ目的変数の値。 | |
X | 適当な関数を説明変数に適用して作ったデザイン行列。確率変数ではない。 | |
X[:, j-1] | デザイン行列から列目をとったもの。 | |
coef | 求めたい回帰係数。推定値だけどハットは省略した。 | |
yhat | 訓練標本に対する推測値。 | |
resid[i-1] | で、推測値と実値の誤差。 | |
eta | 学習の際に用いる閾値。値は適当にセットした。 | |
w_diag[i-1] | 表が崩れてしまったので別途下に記載 | |
len(y) | 訓練標本の数となっている。 | |
W | 対角行列を指す。 | |
- | 行列に対してはの一般逆行列を指す。 |
iris
の呼び出しと、線形回帰係数を得るまで。# インポート
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
coef
を更新する。def makew(ndarr, eta):
w_diag = np.where(ndarr <= eta, ndarr, eta) / ndarr
W = np.diag(w_diag)
return W
def make_abs_resid(y, yhat):
resid = np.abs(y - yhat)
return resid
yhat = X @ coef
resid = make_abs_resid(y, yhat)
eta = 0.5
W = makew(resid, eta)
# 疑似コードに収束するまで、と書いてあるが今回は収束を確認するための閾値を設けずに1万回だけ更新する。
for _ in range(10000):
yhat = X @ coef
resid = make_abs_resid(y, yhat)
W = makew(resid, eta)
coef = np.linalg.inv(X.T @ W @ X) @ X.T @ W @ y
coef