Numpyだけで回帰分析その5。linalgのlstsqで重回帰分析

Numpyだけで回帰分析その5。
どうやら重回帰分析もできるの巻。

Numpyのlinalg.lstsq()、線形代数的に最小2乗法を扱うメソッドを使う。polyfit()・poly1d()と何が違うかは後の課題としておこ。
参考
scipyのofficial document
numpy.linalg.lstsq — NumPy v1.17 Manual

実行環境

---

Androidスマホ
termux
Python3.8
JupyterNotebook

• Pythonライブラリ
  • Numpy
  • Pandas
• その他
  • スマホアプリのExcel

準備するもの

---

ワインの成分の重回帰分析学習用データみたいなものを使うと電話が火を吹いて壊れそうなので、適当にWebページにあった、体重・身長・ウエスト・胸囲の小さなデータをExcelに入力して使うことにする。

4列のデータ。

• 体重 : 目的変数
• その他 : 説明変数

その他から体重を割り出せる回帰モデルができれば成功かな？

PandasでExcelファイルの読み込み

---

import pandas as pd

# Excelファイル
data_file = "body.xlsx"

# PandasでExcelファイルの読み込み
df = pd.read_excel(data_file)

# 表示
df

	体重	身長	腹囲	胸囲
0	82.3	177.5	80.2	95.2
1	70.4	173.0	75.6	83.3
2	65.3	168.8	70.2	78.3
3	62.2	165.5	67.0	75.6
4	54.1	164.5	62.4	72.3
5	50.0	160.2	60.0	69.9
6	58.8	170.2	69.4	77.9
7	77.7	180.2	80.2	88.2
8	92.8	182.3	88.3	93.5
9	60.4	160.3	62.0	74.9

10行4列。これが全て。

# いろんな情報とりあえず
print(type(df)," 【型】\n\n",
      df.shape, " 【形状】\n\n",
      df.size, " 【要素数】\n\n",
      df.info(), " 【info情報】\n\n",
      df.describe(), " 【describe情報】\n\n",
      df.isnull().sum(), " 【欠損値の数】\n\n",
)

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 10 entries, 0 to 9
Data columns (total 4 columns):
体重    10 non-null float64
身長    10 non-null float64
腹囲    10 non-null float64
胸囲    10 non-null float64
dtypes: float64(4)
memory usage: 448.0 bytes
<class 'pandas.core.frame.DataFrame'>  【型】

 (10, 4)  【形状】

 40  【要素数】

 None  【info情報】

               体重          身長         腹囲         胸囲
count  10.000000   10.000000  10.000000  10.000000
mean   67.400000  170.250000  71.530000  80.910000
std    13.401161    7.902637   9.321904   8.800941
min    50.000000  160.200000  60.000000  69.900000
25%    59.200000  164.750000  63.550000  75.075000
50%    63.750000  169.500000  69.800000  78.100000
75%    75.875000  176.375000  79.050000  86.975000
max    92.800000  182.300000  88.300000  95.200000  【describe情報】

 体重    0
身長    0
腹囲    0
胸囲    0
dtype: int64  【欠損値の数】

Numpyのlinalg.lstsqで重回帰分析

---

山場。
説明変数の行列の一番右に「1」で構成する列を追加し、linalg.lstsq()メソッドに渡す。説明変数ごとの係数などが算出される。

【参考】
重回帰分析をするための何種類かの方法が書いてあって参考になりました。

https://www.ydc.co.jp/column/0002/20181207l.html

import numpy as np

# 体重 目的変数  二次元に
y = df.iloc[:,0].values.reshape(-1,1)

# 身長ウエスト胸囲 説明変数 
x = df.values[:,1:]

# 型とオブジェクトの表示(省略)
#print(type(y),type(x))
print(y)
print("")
print(x)

[[82.3]
 [70.4]
 [65.3]
 [62.2]
 [54.1]
 [50. ]
 [58.8]
 [77.7]
 [92.8]
 [60.4]]

[[177.5  80.2  95.2]
 [173.   75.6  83.3]
 [168.8  70.2  78.3]
 [165.5  67.   75.6]
 [164.5  62.4  72.3]
 [160.2  60.   69.9]
 [170.2  69.4  77.9]
 [180.2  80.2  88.2]
 [182.3  88.3  93.5]
 [160.3  62.   74.9]]

# 説明変数の行数
n = x.shape[0] 
# 行数表示
n

10

# 定数項 二次元の1だけの列ベクトル生成
ones = np.ones(n).reshape(-1,1)
print(ones)

[[1.]
 [1.]
 [1.]
 [1.]
 [1.]
 [1.]
 [1.]
 [1.]
 [1.]
 [1.]]

# 定数項を横に結合
x = np.hstack([x, ones])
# 形状表示
print(x.shape)
print(x)

(10, 5)
[[177.5  80.2  95.2   1.   ]
 [173.   75.6  83.3   1.    ]
 [168.8  70.2  78.3   1.   ]
 [165.5  67.   75.6   1.    ]
 [164.5  62.4  72.3   1.   ]
 [160.2  60.   69.9   1.    ]
 [170.2  69.4  77.9   1.   ]
 [180.2  80.2  88.2   1.   ]
 [182.3  88.3  93.5   1.   ]
 [160.3  62.   74.9   1.   ]]

# 回帰多項式の係数生成
coef = np.linalg.lstsq(x, y ,rcond=None)

# 表示
print(coef)

(array([[ -1.28680201],
       [  1.90992916],
       [  0.61255179],
       [100.29924438]]), array([32.92019628]), 4, array([6.38817982e+02, 2.01529425e+01, 6.08380411e+00, 4.28781259e-02]))

注意
linalg.lstsq()の引数に「rcond=None」と書いています。 警告が出たので書きました。
小さなシンギュラーバリューのカットオフ比率の設定云々と。意味が理解できないが「警告を黙らせたいならこうしろ、古いやつ使うならrcond=-1を使え」とは書いてある。
[https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html:title]

回帰係数算出その後

---

一応ここまでが、Numpyのlinalg.lstsq()で重回帰分析の値を算出するまでの手順は完了です。

y = a1x1 + a2x2 + a3x3 + …

みたいな多項式の係数aが配列で算出されました。
あとは戻り値が何なのかの確認と、実際の値を再代入し体重が予測出来ているかを我田引水、自作自演、自家薬籠中的ではありますが確認します。

linalg.lstsqの戻り値

---

スーパースター、モドリッチ。

さて、変数名「coef」として取り出した値だが、傾き以外にも何種類か返されている。

• 最小2乗解の配列。多項式の係数の配列
• residuals 残差の合計の配列
• 行列ランクを整数型
• 係数の特異値
  
  
• 警告
  • LinAlgError 計算が収束しない場合に発生

予測モデルに再代入

---

係数を出すために使った説明変数x、身長ウエスト胸囲を多項式に代入して体重が出せるかを確認する。

print("・ナンバー0某氏のx\n 身長/ウエスト/胸囲\n",  x[0])
print(" 体重y\n", y[0])
print("\n・係数の個数", coef[0].size)
print("・係数\n", coef[0])

・ナンバー0某氏のx
 身長/ウエスト/胸囲
 [177.5  80.2  95.2   1.    1. ]
 体重y
 [82.3]

・係数の個数 4
・係数
 [[ -1.28680201]
 [  1.90992916]
 [  0.61255179]
 [100.29924438]]

print("0氏の体重 :",y[0])

# 0氏の説明変数
x0 = x[0]

# 体重予測モデルに0氏を当てはめ
y = coef[0][0]*x0[0] + coef[0][1]*x0[1] + coef[0][2]*x0[2]  + coef[0][3]
print("0氏の予測体重 ",y)

0氏の体重 : [82.3]
0氏の予測体重  [83.38313627]

まぁまぁ合っている。
未知の人物の身長ウエスト胸囲のデータをモデルに代入すれば、その人の体重が大体予測できるという事でよろしいかな。

おわりに

---

意外と簡単に、Numpyだけで重回帰分析ができる事が分かった。
定数項の1だけの列を入れる理由や仕組みは後程調べるとして、問題はこれで得られる係数がsklearnやtensorflowなどと比較して精度がどれ程違うか、ですかね。

ほぼ写経だったが、やる価値のあるものだった。
次はスマホが火を吹きそうなデータで試してみようww。