Ccmmutty logo
Commutty IT
・4 min read

スパース学習についてその1 L1制約付き最小二乗学習

Python機械学習スパース推定
https://cdn.magicode.io/media/notebox/73fe7985-eaa3-427e-a9d8-1493b0edfe7a.jpeg

参考文献

文献1.のChapter5、特に5.2を参考としている。
  1. https://bookclub.kodansha.co.jp/product?item=0000148211

L1制約付き最小二乗学習

概要

説明変数により構成されたデザイン行列Φ\Phi、目的変数yyについてスパース学習を行う方法の1つであるl1l_1制約付き最小二乗学習について紹介する。 その方法は、推定するパラメータθ\thetaのL1ノルム∣∣θ∣∣1\mid\mid \theta \mid\mid_1に対して制約条件をつけた上でθ⊤\theta^{\top}と二乗誤差∣∣y−Φθ∣∣22\mid\mid y-\Phi\theta \mid\mid^2_2を最小化することである。表式にすると、以下のようになる。
JLS(θ)=∣∣y−Φθ∣∣220<R,R∈Rmin⁡θJLS(θ)s.t.∣∣θ∣∣1≤RJ_{LS}\left( \theta \right) = \mid\mid y-\Phi\theta \mid\mid^2_2 \\ 0 < R, R \in \mathbb{R} \\ \min_{\theta} J_{LS}\left( \theta \right) \text{s.t.} \mid\mid \theta \mid\mid_1 \leq R
このように表式にはしたが、これをこのまま最小化するのではなく等価で解きやすい別の問題に置き換えて解くことによってパラメータを推定する方法について以下で紹介する。

解く問題について

最終的には、絶対値関数∣θ∣\mid \theta \midが微分できる二次関数θ22c+c2\frac{\theta^2}{2c}+\frac{c}{2}で上からおさえられる事実を使って、以下の問題を解くことになる。この問題は、λR\lambda_RとRRについて整理すれば概要の部分で表式とした問題と同じものとなっている。解き方については、疑似コードの部分で計算方法を紹介する。
JLS(θ)=∣∣y−Φθ∣∣22J=JLS(θ)+λR∣∣θ∣∣1J_{LS}\left( \theta \right) = \mid\mid y-\Phi\theta \mid\mid^2_2 \\ J = J_{LS}\left( \theta \right) + \lambda_R\mid\mid \theta \mid\mid_1 \\

記号とパラメータ

数式中の記号コード中の変数概要
yyy目的変数で、実数値をとる。
yiy_iy[i-1]ii番目の標本が持つ目的変数の値。
Φ\PhiX関数Φ\Phiを説明変数XXに適用したデザイン行列。
ϕ(xj)\phi(x_j)X[:, j-1]デザイン行列からjj列目をとったもの。
θ\thetacoef求めたい回帰係数。推定値。
fθ(xi)f_{\theta}(x_i)yhat訓練標本iiに対する推測値。
rir_iresid[i-1]yi−fθ(xi)y_i - f_{\theta}(x_i)で、推測値と実値の誤差。
λR\lambda_Rlmd制約条件のRRによって決定されるハイパーパラメータ
θi\theta_icoef[i-1]推定したii番目のパラメータθi\theta_i
bblen(coef)推定する回帰係数のパラメータ数
Θ\ThetaTheta対角行列diag(∣θ1∣,…,∣θb∣)diag(\mid \theta_1 \mid,\dots,\mid \theta_b \mid)を指す。
†\dagger-行列AAに対してA†A^{\dagger}はAAの一般逆行列を指す。

疑似コード

  1. θ\thetaの初期値を適当に定める。本記事では(Φ⊤Φ)†Φ⊤y\left( \Phi^{\top}\Phi \right)^{\dagger}\Phi^{\top}yとした。つまり、普通に線形回帰係数を得た。
  2. 推定値θ\thetaの値が収束するまで、以下の2つのステップを繰り替えす。本記事では、値の収束ではなく10000回計算することとした。
    1. 現在の推定値θ\thetaを使って、対角行列Θ\Thetaを計算する。
    2. θ\thetaを(Φ⊤Φ+λRΘ†)−1Φ⊤y\left( \Phi^{\top}\Phi + \lambda_R \Theta^{\dagger} \right)^{-1}\Phi^{\top}yで更新する。

データセットirisの呼び出しと、線形回帰係数を得るまで。

python
!pip install statsmodels
Collecting statsmodels
Downloading statsmodels-0.13.2-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (9.8 MB) [?25l | | 10 kB 15.5 MB/s eta 0:00:01 | | 20 kB 18.2 MB/s eta 0:00:01 | | 30 kB 10.1 MB/s eta 0:00:01 |▏ | 40 kB 8.4 MB/s eta 0:00:02 |▏ | 51 kB 6.9 MB/s eta 0:00:02 |▏ | 61 kB 8.1 MB/s eta 0:00:02 |▎ | 71 kB 8.2 MB/s eta 0:00:02 |▎ | 81 kB 9.1 MB/s eta 0:00:02 |▎ | 92 kB 8.3 MB/s eta 0:00:02 |▎ | 102 kB 8.0 MB/s eta 0:00:02 |▍ | 112 kB 8.0 MB/s eta 0:00:02 |▍ | 122 kB 8.0 MB/s eta 0:00:02 |▍ | 133 kB 8.0 MB/s eta 0:00:02 |▌ | 143 kB 8.0 MB/s eta 0:00:02 |▌ | 153 kB 8.0 MB/s eta 0:00:02 |▌ | 163 kB 8.0 MB/s eta 0:00:02 |▋ | 174 kB 8.0 MB/s eta 0:00:02 |▋ | 184 kB 8.0 MB/s eta 0:00:02 |▋ | 194 kB 8.0 MB/s eta 0:00:02 |▋ | 204 kB 8.0 MB/s eta 0:00:02 |▊ | 215 kB 8.0 MB/s eta 0:00:02 |▊ | 225 kB 8.0 MB/s eta 0:00:02 |▊ | 235 kB 8.0 MB/s eta 0:00:02 |▉ | 245 kB 8.0 MB/s eta 0:00:02 |▉ | 256 kB 8.0 MB/s eta 0:00:02 |▉ | 266 kB 8.0 MB/s eta 0:00:02 |█ | 276 kB 8.0 MB/s eta 0:00:02 |█ | 286 kB 8.0 MB/s eta 0:00:02 |█ | 296 kB 8.0 MB/s eta 0:00:02 |█ | 307 kB 8.0 MB/s eta 0:00:02 |█ | 317 kB 8.0 MB/s eta 0:00:02 |█ | 327 kB 8.0 MB/s eta 0:00:02 |█ | 337 kB 8.0 MB/s eta 0:00:02 |█▏ | 348 kB 8.0 MB/s eta 0:00:02 |█▏ | 358 kB 8.0 MB/s eta 0:00:02 |█▏ | 368 kB 8.0 MB/s eta 0:00:02 |█▎ | 378 kB 8.0 MB/s eta 0:00:02 |█▎ | 389 kB 8.0 MB/s eta 0:00:02 |█▎ | 399 kB 8.0 MB/s eta 0:00:02 |█▎ | 409 kB 8.0 MB/s eta 0:00:02 |█▍ | 419 kB 8.0 MB/s eta 0:00:02 |█▍ | 430 kB 8.0 MB/s eta 0:00:02 |█▍ | 440 kB 8.0 MB/s eta 0:00:02 |█▌ | 450 kB 8.0 MB/s eta 0:00:02 |█▌ | 460 kB 8.0 MB/s eta 0:00:02 |█▌ | 471 kB 8.0 MB/s eta 0:00:02 |█▋ | 481 kB 8.0 MB/s eta 0:00:02 |█▋ | 491 kB 8.0 MB/s eta 0:00:02 |█▋ | 501 kB 8.0 MB/s eta 0:00:02 |█▋ | 512 kB 8.0 MB/s eta 0:00:02 |█▊ | 522 kB 8.0 MB/s eta 0:00:02 |█▊ | 532 kB 8.0 MB/s eta 0:00:02 |█▊ | 542 kB 8.0 MB/s eta 0:00:02 |█▉ | 552 kB 8.0 MB/s eta 0:00:02 |█▉ | 563 kB 8.0 MB/s eta 0:00:02 |█▉ | 573 kB 8.0 MB/s eta 0:00:02 |██ | 583 kB 8.0 MB/s eta 0:00:02 |██ | 593 kB 8.0 MB/s eta 0:00:02 |██ | 604 kB 8.0 MB/s eta 0:00:02 |██ | 614 kB 8.0 MB/s eta 0:00:02 |██ | 624 kB 8.0 MB/s eta 0:00:02 |██ | 634 kB 8.0 MB/s eta 0:00:02 |██ | 645 kB 8.0 MB/s eta 0:00:02 |██▏ | 655 kB 8.0 MB/s eta 0:00:02 |██▏ | 665 kB 8.0 MB/s eta 0:00:02 |██▏ | 675 kB 8.0 MB/s eta 0:00:02 |██▎ | 686 kB 8.0 MB/s eta 0:00:02 |██▎ | 696 kB 8.0 MB/s eta 0:00:02 |██▎ | 706 kB 8.0 MB/s eta 0:00:02 |██▎ | 716 kB 8.0 MB/s eta 0:00:02
|██▍ | 727 kB 8.0 MB/s eta 0:00:02 |██▍ | 737 kB 8.0 MB/s eta 0:00:02 |██▍ | 747 kB 8.0 MB/s eta 0:00:02 |██▌ | 757 kB 8.0 MB/s eta 0:00:02 |██▌ | 768 kB 8.0 MB/s eta 0:00:02 |██▌ | 778 kB 8.0 MB/s eta 0:00:02 |██▋ | 788 kB 8.0 MB/s eta 0:00:02 |██▋ | 798 kB 8.0 MB/s eta 0:00:02 |██▋ | 808 kB 8.0 MB/s eta 0:00:02 |██▋ | 819 kB 8.0 MB/s eta 0:00:02 |██▊ | 829 kB 8.0 MB/s eta 0:00:02 |██▊ | 839 kB 8.0 MB/s eta 0:00:02 |██▊ | 849 kB 8.0 MB/s eta 0:00:02 |██▉ | 860 kB 8.0 MB/s eta 0:00:02 |██▉ | 870 kB 8.0 MB/s eta 0:00:02 |██▉ | 880 kB 8.0 MB/s eta 0:00:02 |███ | 890 kB 8.0 MB/s eta 0:00:02 |███ | 901 kB 8.0 MB/s eta 0:00:02 |███ | 911 kB 8.0 MB/s eta 0:00:02 |███ | 921 kB 8.0 MB/s eta 0:00:02 |███ | 931 kB 8.0 MB/s eta 0:00:02 |███ | 942 kB 8.0 MB/s eta 0:00:02 |███ | 952 kB 8.0 MB/s eta 0:00:02 |███▏ | 962 kB 8.0 MB/s eta 0:00:02 |███▏ | 972 kB 8.0 MB/s eta 0:00:02 |███▏ | 983 kB 8.0 MB/s eta 0:00:02 |███▎ | 993 kB 8.0 MB/s eta 0:00:02 |███▎ | 1.0 MB 8.0 MB/s eta 0:00:02 |███▎ | 1.0 MB 8.0 MB/s eta 0:00:02 |███▎ | 1.0 MB 8.0 MB/s eta 0:00:02 |███▍ | 1.0 MB 8.0 MB/s eta 0:00:02 |███▍ | 1.0 MB 8.0 MB/s eta 0:00:02 |███▍ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▌ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▌ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▌ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▋ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▋ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▋ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▋ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▊ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▊ | 1.1 MB 8.0 MB/s eta 0:00:02 |███▊ | 1.2 MB 8.0 MB/s eta 0:00:02 |███▉ | 1.2 MB 8.0 MB/s eta 0:00:02 |███▉ | 1.2 MB 8.0 MB/s eta 0:00:02 |███▉ | 1.2 MB 8.0 MB/s eta 0:00:02 |████ | 1.2 MB 8.0 MB/s eta 0:00:02 |████ | 1.2 MB 8.0 MB/s eta 0:00:02 |████ | 1.2 MB 8.0 MB/s eta 0:00:02 |████ | 1.2 MB 8.0 MB/s eta 0:00:02 |████ | 1.2 MB 8.0 MB/s eta 0:00:02 |████ | 1.2 MB 8.0 MB/s eta 0:00:02 |████ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▏ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▏ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▏ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▎ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▎ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▎ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▎ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▍ | 1.3 MB 8.0 MB/s eta 0:00:02 |████▍ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▍ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▌ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▌ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▌ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▋ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▋ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▋ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▋ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▊ | 1.4 MB 8.0 MB/s eta 0:00:02 |████▊ | 1.5 MB 8.0 MB/s eta 0:00:02 |████▊ | 1.5 MB 8.0 MB/s eta 0:00:02 |████▉ | 1.5 MB 8.0 MB/s eta 0:00:02 |████▉ | 1.5 MB 8.0 MB/s eta 0:00:02 |████▉ | 1.5 MB 8.0 MB/s eta 0:00:02 |█████ | 1.5 MB 8.0 MB/s eta 0:00:02 |█████ | 1.5 MB 8.0 MB/s eta 0:00:02 |█████ | 1.5 MB 8.0 MB/s eta 0:00:02 |█████ | 1.5 MB 8.0 MB/s eta 0:00:02 |█████ | 1.5 MB 8.0 MB/s eta 0:00:02 |█████ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▏ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▏ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▏ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▎ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▎ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▎ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▎ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▍ | 1.6 MB 8.0 MB/s eta 0:00:02 |█████▍ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▍ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▌ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▌ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▌ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▋ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▋ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▋ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▋ | 1.7 MB 8.0 MB/s eta 0:00:02 |█████▊ | 1.8 MB 8.0 MB/s eta 0:00:02 |█████▊ | 1.8 MB 8.0 MB/s eta 0:00:02 |█████▊ | 1.8 MB 8.0 MB/s eta 0:00:02 |█████▉ | 1.8 MB 8.0 MB/s eta 0:00:02 |█████▉ | 1.8 MB 8.0 MB/s eta 0:00:02 |█████▉ | 1.8 MB 8.0 MB/s eta 0:00:02 |█████▉ | 1.8 MB 8.0 MB/s eta 0:00:02 |██████ | 1.8 MB 8.0 MB/s eta 0:00:02 |██████ | 1.8 MB 8.0 MB/s eta 0:00:02 |██████ | 1.8 MB 8.0 MB/s eta 0:00:02 |██████ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████▏ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████▏ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████▏ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████▏ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████▎ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████▎ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████▎ | 1.9 MB 8.0 MB/s eta 0:00:01 |██████▍ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▍ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▍ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▌ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▌ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▌ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▌ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▋ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▋ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▋ | 2.0 MB 8.0 MB/s eta 0:00:01 |██████▊ | 2.1 MB 8.0 MB/s eta 0:00:01 |██████▊ | 2.1 MB 8.0 MB/s eta 0:00:01 |██████▊ | 2.1 MB 8.0 MB/s eta 0:00:01 |██████▉ | 2.1 MB 8.0 MB/s eta 0:00:01 |██████▉ | 2.1 MB 8.0 MB/s eta 0:00:01 |██████▉ | 2.1 MB 8.0 MB/s eta 0:00:01 |██████▉ | 2.1 MB 8.0 MB/s eta 0:00:01 |███████ | 2.1 MB 8.0 MB/s eta 0:00:01 |███████ | 2.1 MB 8.0 MB/s eta 0:00:01 |███████ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████▏ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████▏ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████▏ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████▏ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████▎ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████▎ | 2.2 MB 8.0 MB/s eta 0:00:01 |███████▎ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▍ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▍ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▍ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▌ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▌ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▌ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▌ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▋ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▋ | 2.3 MB 8.0 MB/s eta 0:00:01 |███████▋ | 2.4 MB 8.0 MB/s eta 0:00:01 |███████▊ | 2.4 MB 8.0 MB/s eta 0:00:01 |███████▊ | 2.4 MB 8.0 MB/s eta 0:00:01 |███████▊ | 2.4 MB 8.0 MB/s eta 0:00:01 |███████▉ | 2.4 MB 8.0 MB/s eta 0:00:01 |███████▉ | 2.4 MB 8.0 MB/s eta 0:00:01 |███████▉ | 2.4 MB 8.0 MB/s eta 0:00:01 |███████▉ | 2.4 MB 8.0 MB/s eta 0:00:01 |████████ | 2.4 MB 8.0 MB/s eta 0:00:01 |████████ | 2.4 MB 8.0 MB/s eta 0:00:01 |████████ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████▏ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████▏ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████▏ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████▏ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████▎ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████▎ | 2.5 MB 8.0 MB/s eta 0:00:01 |████████▎ | 2.6 MB 8.0 MB/s eta 0:00:01 |████████▍ | 2.6 MB 8.0 MB/s eta 0:00:01
|████████▍ | 2.6 MB 8.0 MB/s eta 0:00:01 |████████▍ | 2.6 MB 8.0 MB/s eta 0:00:01 |████████▌ | 2.6 MB 8.0 MB/s eta 0:00:01 |████████▌ | 2.6 MB 8.0 MB/s eta 0:00:01 |████████▌ | 2.6 MB 8.0 MB/s eta 0:00:01 |████████▌ | 2.6 MB 8.0 MB/s eta 0:00:01 |████████▋ | 2.6 MB 8.0 MB/s eta 0:00:01 |████████▋ | 2.7 MB 8.0 MB/s eta 0:00:01 |████████▋ | 2.7 MB 8.0 MB/s eta 0:00:01 |████████▊ | 2.7 MB 8.0 MB/s eta 0:00:01 |████████▊ | 2.7 MB 8.0 MB/s eta 0:00:01 |████████▊ | 2.7 MB 8.0 MB/s eta 0:00:01 |████████▉ | 2.7 MB 8.0 MB/s eta 0:00:01 |████████▉ | 2.7 MB 8.0 MB/s eta 0:00:01 |████████▉ | 2.7 MB 8.0 MB/s eta 0:00:01 |████████▉ | 2.7 MB 8.0 MB/s eta 0:00:01 |█████████ | 2.7 MB 8.0 MB/s eta 0:00:01 |█████████ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████▏ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████▏ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████▏ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████▏ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████▎ | 2.8 MB 8.0 MB/s eta 0:00:01 |█████████▎ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▎ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▍ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▍ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▍ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▌ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▌ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▌ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▌ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▋ | 2.9 MB 8.0 MB/s eta 0:00:01 |█████████▋ | 3.0 MB 8.0 MB/s eta 0:00:01 |█████████▋ | 3.0 MB 8.0 MB/s eta 0:00:01 |█████████▊ | 3.0 MB 8.0 MB/s eta 0:00:01 |█████████▊ | 3.0 MB 8.0 MB/s eta 0:00:01 |█████ too many strings
python
# インポート
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])

実際にcoefを更新する。

python
def make_Theta(coef):
  abs_coef = np.abs(coef)
  Theta = np.diag(abs_coef)
  return Theta
python
lmd = 10
python
# 疑似コードに収束するまで、と書いてあるが今回は収束を確認するための閾値を設けずに1万回だけ更新する。
python
for _ in range(10000):
    Theta = make_Theta(coef)
    coef = np.linalg.inv(X.T @ X + lmd * np.linalg.pinv(Theta)) @ X.T @ y
python
coef
array([ 4.54049476e-05, 7.46181973e-01, -3.54373634e-01])

Discussion

コメントにはログインが必要です。