Python初心者の備忘録 #15 ~機械学習入門編01~

# はじめに
今回私は最近はやりのchatGPTに興味を持ち、深層学習について学んでみたいと思い立ちました！  
深層学習といえばPythonということなので、最終的にはPythonを使って深層学習ができるとこまでコツコツと学習していくことにしました。  
ただ、勉強するだけではなく少しでもアウトプットをしようということで、備忘録として学習した内容をまとめていこうと思います。  
この記事が少しでも誰かの糧になることを願っております！
**※投稿主の環境はWindowsなのでMacの方は多少違う部分が出てくると思いますが、ご了承ください。**
最初の記事：[Python初心者の備忘録 #01](https://qiita.com/Yushin-Tati/items/961dc47d6163f944f7e9)
前の記事：[Python初心者の備忘録 #14 ~統計学入門編04~](https://qiita.com/Yushin-Tati/items/8aeb2c190e3494dd6b78)
次の記事：[Python初心者の備忘録 #16 ~機械学習入門編02~](https://qiita.com/Yushin-Tati/items/5fe645f0a3c70b4dfa1a)

本記事からは機械学習の入門的なことについてまとめていこうと思います。
今回は**線形回帰**や**その解釈**、**正規方程式**、**特徴量スケーリング**、**ダミー変数**についてまとめております。

:::note warn
本記事はDockerで環境構築を行っているので、環境構築がまだという方は下記記事を参考にして、環境構築をしてください。
[Python初心者の備忘録 #06 ~DSに使われるライブラリ編01~](https://qiita.com/Yushin-Tati/items/5dcad7518a0f8a3fdbc1)

上記に加えて、下記コマンドを実行してpullしてきたimage名にdocker-compose.ymlを修正してください。
```
docker pull datascientistus/ds-python-env3
```

※また2024/04/13以前に環境構築された方は、Dockerfileの内容を更新しているので自身の環境の設定ファイルの更新もお願いいたします。
:::
## ■学習に使用している資料
Udemy：[【前編】米国データサイエンティストがやさしく教える機械学習超入門【Pythonで実践】
](https://www.udemy.com/course/mlpython-1/)
## ■線形回帰(linear regression)
:::note info
今回扱うデータは基本的に標本で、**背後に母集団がある**ということを頭の片隅に入れておいてください！
:::
### ▶線形回帰とは
- 機械学習には大きく分けて回帰と分類(クラスタリング)というのがあり、回帰はもともとあるデータから推測するという手法で、線形回帰はその一つ

例 家賃と広さの関係
![ex_linear_regression.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/ea888918-1c9f-2a05-2736-5d4ac9552a27.jpeg)
上図であれば家賃と広さの関係は次のような式で表すことができる。
$Y=f(X_1)+\epsilon$
\- $Y$：真の目的変数(ex.家賃)の値
\- $X_1$：特徴量(ex.広さ) <- 複数になることもある
\- $f$：真の関数 <- 特徴量と目的変数の関係を表す式
\- $\epsilon$：誤差項 <- $X_1$だけでは説明できない誤差(ex.築年数、駅からの距離、設備)
- 本来であれは上記の式で求めたいがわからないので、機械学習(線形回帰)で真の関数$f$を推測する($\hat{f}$を求める)
$\hat{Y}=\hat{f}(X)$
\- $\hat{Y}$：予測値
\- $\hat{f}$：真の関数の推測値
- 線形回帰では$f$が「線形」と仮定して「**最小二乗法**」で最適なパラメータ($\theta$)を見つける
$f(X)=\theta_0+\theta_1X_1+...+\theta_nX_n$ ($\theta_n$：n番目のパラメータ、$X_n$：n番目の特徴量の値)
:::note info
機械学習では最適なパラメータを推測することを「**学習する**」という。
:::
### ▶最小二乗法(least squares)
- ある特徴量と目的変数の関係を表す関数が直線だと仮定した際に、パラメータ$\theta$を求めるときに採用される方法
- 直線(関数)とそれぞれの標本との誤差(残差)の二乗和が最小になるような$\theta$を求める

$\sum_{i=1}^me_i^2=\sum_{i=1}^m\\{y_i-(\theta_0+\theta_1x_i)\\}^2$
($y_i$：i番目の実データ、m：データの数、$x_i$：ある特徴量のi番目の値)
- 上記の式における$\theta_0$と$\theta_1$を**最急降下法**によって求める

![explain_least_squares.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/7d050d8a-09cf-24f0-ac06-b87085ff1340.jpeg)

### ▶最急降下法(gradient descent)
- 関数の最小値を探索するアルゴリズムの一つ
- 機械学習では最小化する対象となる関数のことを「損失関数(Loss function)」という
※残差のことを損失と捉えている

**残差の二乗の平均(Mean Squared Error：MSE)**
$L(\theta_0, \theta_1)=\frac{1}{m}\sum_{i=1}^me_i^2=\frac{1}{m}\sum_{i=1}^m\\{y_i-(\theta_0+\theta_1x_i)\\}^2$
![MSE_3D_graph.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/52186bd8-d854-aecb-a8b4-9a6a2e8d7adb.jpeg)

- 関数が最小になる変数の値を探索することを最適化問題といい、最急降下法はその中の勾配法に分類される
- 変数$\theta$における関数の接線の傾きを勾配ベクトルといい、勾配ベクトルが0になる方向に少しずつ$\theta$を調整する
$\theta:=\theta-\alpha\frac{d}{d\theta}L(\theta)$ ($\alpha$：学習率)
※学習率$\alpha$が大きいほど探索が急激に進み、小さいほど緩やかに進む

![example_image_gradient_descent.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/8395c2d5-2803-f951-a226-397dbb433e54.jpeg)
![example_image_gradient_descent_2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/046781c4-d267-c2dc-e3f9-f81039bdf878.jpeg)
- $\theta$が複数あったとしてもそれぞれの$\theta$について上記の式を適用し、$\theta$の値を**同時に**更新していく
$\theta_0:=\theta_0-\alpha\frac{\partial}{\partial\theta_0}L(\theta_0,\theta_1)$
$\theta_1:=\theta_1-\alpha\frac{\partial}{\partial\theta_1}L(\theta_0,\theta_1)$
$\frac{\partial}{\partial\theta_0}L(\theta_0,\theta_1) = \frac{2}{m}\sum_{i=1}^m(\theta_0+\theta_1x_i-y_i)$
$\frac{\partial}{\partial\theta_1}L(\theta_0,\theta_1) = \frac{2}{m}\sum_{i=1}^m(\theta_0+\theta_1x_i-y_i)x_i$
※$\frac{\partial}{\partial x}$は偏微分といわれるものだが、複数パラメータに対する微分と考えておけばOK

上記の式にMSEの式を代入すると...(合成関数の微分を使用)
※合成関数の微分：$\\{f(g(x))\\}'=f(g(x))'・g(x)'$
-> 今回であれば$g(x)=y_i-(\theta_0+\theta_1x\_i)$となる
### ▶Pythonで最急降下法
:::note info
実際に手を動かしながら記事を読みたい場合は、下記githubからsample_data.csvをDLしてください。
https://github.com/Yushin-Tati/Learning_machine_learning
:::
#### データの準備
```python
import pandas as pd

df = pd.read_csv('data/sample_data.csv')
```
#### データの可視化
```python
import seaborn as sns

sns.scatterplot(x='space', y='rent', data=df)
```
![scatter_plot_gradient_descent.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/953b1d64-4e63-3355-2935-1ff369d66390.jpeg)
#### 損失関数を計算する関数の定義
```python
import nampy as np

x = df['space'].values
y = df{'rent'}.values

# 損失関数の結果を返す関数
def cost_func(theta_0, theta_1, x, y):
    return np.mean(np.square(y - (theta_0 + theta_1 * x)))
cost_func(5, 7, x=x, y=y)  # -> 138703.20244638805
```
#### 損失関数の可視化
```python
import matplotlib.pyplot as plt
from itertools import product

# 軸の値の数
n0 = n1 = 300

# それぞれのパラメータθの最小値と最大値
min_0 = -20
max_0 = 20
min_1 = -1
max_1 = 1

# 軸の値設定
theta_0 = np.linspace(min_0, max_0, n0)
theta_1 = np.linspace(min_1, max_1, n1)

# plot_surface用に2次元にする
theta_0_ax, theta_1_ax = np.meshgrid(theta_0, theta_1)

# 損失関数の取得
z = [cost_func(param[0], param[1], x=x, y=y) for param in list(product(theta_0_ax, theta_1_ax))]
# plot_surface用に2次元にする
Z = np.array(z).reshape(n0, n1)

# 3Dのグラフになるようにする
ax = plt.axes(projection='3d')
# meshgridとproductおよびreshapeの処理がそれぞれ行と列が逆になっているので転置.Tする
ax.plot_surface(theta_0_ax.T, theta_1_ax.T, Z, cmap='jet', alpha=0.5)

# グラフの見え方(角度やどこから見るか、軸の位置など)を調整
plt.gca().invert_xaxis()
ax.view_init(elev=30, azim=10)
```
![3d_gradient_descent.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/c3b35a19-3a40-003c-8c20-624e85810e7d.jpeg)
#### 最急降下法
1. $\theta_0$と$\theta_1$の初期値を決め、$\theta$更新用の関数を定義する
1. イテレーション数と学習率$\alpha$を設定し、for文で回して$\theta$を更新する
※イテレーション：10万回、$\alpha$：0.00005
1. イテレーションの履歴は保存する(後で推移を確認するため)

```python
# 本来であれば初期値はランダム
# theta_0_init = np.random.uniform(min_0, max_0)
# theta_1_init = np.random.uniform(min_1, max_1)

# 描画時にわかりやすいところに初期値を設定
theta_0_init = -5
theta_1_init = -0.5

# パラメータ更新用の関数を定義
def update_theta0(theta_0, theta_1, x, y, alpha=0.05):
    return theta_0 - alpha * 2 * np.mean((theta_0 + theta_1 * x) - y)

def update_theta1(theta_0, theta_1, x, y, alpha=0.05):
    return theta_1 - alpha * 2 * np.mean(((theta_0 + theta_1 * x) - y) * x) 

# 機械学習では学習の単位はエポック(epoch)
epochs = 100000
# 学習率
alpha = 0.00005

# パラメータ更新の履歴
theta_0_hist = []
theta_1_hist = []
# 初期化
theta_0_hist.append(theta_0_init)
theta_1_hist.append(theta_1_init)

# 最急降下法
for _ in range(epochs):
    updated_theta_0 = update_theta0(theta_0_hist[-1], theta_1_hist[-1], x=x, y=y, alpha=alpha)
    updated_theta_1 = update_theta1(theta_0_hist[-1], theta_1_hist[-1], x=x, y=y, alpha=alpha)
    theta_0_hist.append(updated_theta_0)
    theta_1_hist.append(updated_theta_1)

# 学習結果
print(theta_0_hist[-1])  # -> 5.692313583745944
print(theta_1_hist[-1])  # -> 0.2972194648547114
```
#### 結果を描画
1. それぞれの$\theta_0$と$\theta_1$における損失関数の推移を計算し、損失関数に合わせて描画
1. $\theta_0$と$\theta_1$の推移を等高線で描画
1. 得られた$\theta_0$と$\theta_1$の最適解で回帰直線を描画

```python
# 損失の推移
cost_hist = [cost_func(*param, x=x, y=y) for param in zip(theta_0_hist, theta_1_hist)]

# 3D描画
ax = plt.axes(projection='3d')
ax.plot(theta_0_hist, theta_1_hist, cost_hist, 'x-')
ax.plot_surface(theta_0_ax.T, theta_1_ax.T, Z, cmap='jet', alpha=0.5)
plt.gca().invert_xaxis()
ax.view_init(elev=30, azim=10)
```
![3d_transition_gradient_descent.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/1d2fa92e-435c-6694-7130-ad88efa2745d.jpeg)
```python
# 2D描画(等高線)
plt.contour(theta_1_ax.T, theta_0_ax.T, Z, 100, cmap='jet')
plt.plot(theta_1_hist, theta_0_hist, 'x-')
plt.xlabel('theta 1')
plt.ylabel('theta 0')
```
![2d_transition_gradient_descent.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/c798b152-737b-3d48-cb4c-540867cf8874.jpeg)
```python
# 回帰直線
sns.scatterplot(x='space', y='rent', data=df)
x_values = np.arange(120)
y_values = theta_0_hist[-1] + theta_1_hist[-1]*x_values
plt.plot(x_values, y_values, '-', color='r')
```
![result_linear_regression.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/9675fe36-ee7e-95e2-2b27-3cf3d19309b1.jpeg)
```python
# 回帰直線を使って予測値を求める
space = 70
rent = theta_0_hist[-1] + theta_1_hist[-1]*space
rent  # -> 26.497676123575744
```
### ▶最急降下法の注意点
- 設定した初期値によっては、最適解(global optima)ではなく、局所解(local optima)にたどり着くことがある
⇒何回か違う初期値を設定して、最も値が低くなる解を導き出すようにする
- 学習率$\alpha$が大きすぎると、発散して学習が進まない
⇒$\alpha$は重要なハイパーパラメータなので、適切な値を設定する手腕を磨く

![warning_gradient_descent.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/a4ce7dd8-00fe-2dd2-a9b6-54d68a7376de.jpeg)
## ■正規方程式(normal equation)
#### ※行列計算
簡単な復習程度で...
![matrix_calculation_1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/e3314aba-b5cf-8c7c-c576-32b2fc0de044.jpeg)
![matrix_calculation_2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/0f40108f-1aa2-9e4e-7fae-1c29faebec8d.jpeg)
### ▶正規方程式とは
- 最急降下法を使用せずに線形回帰のパラメータ$\theta$を解析的の求める方程式で、線形回帰のモデルでのみ使用することができる
- 一発で確実に買いを求めることが可能

損失関数の式を行列で表すと下図のように変換することができる。
![lost_change_to_matrix_calculation.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/bf30e4d3-66df-a378-21a1-da5a9ca86026.jpeg)
※$\hat{y}=\theta_0+\theta_1x_{i1}$を行列で表すと、$\hat{y}=\begin{bmatrix}1 & x_{i1}\end{bmatrix}\begin{bmatrix}\theta_0 \\\\ \theta_1 \end{bmatrix}$となることから$\hat{y}$は上記のようにあらわすことができる。

よって、求めたい損失関数が最小となるパラメータ$\theta$の最適解は、下図のようになる。
![normal_equation.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/77a5f780-8a36-f49d-b6e2-f959561a9d95.jpeg)
### ▶Pythonで正規方程式
1. $X$を作成
1. $y$と$X$を使って、正規方程式で$\theta$を求める
$\theta=(X^TX)^{-1}X^Ty$
1. 結果を最急降下法と比較する形で等高線で描画
#### Xを作成
```python
X = np.vstack([np.ones(len(x)), x]).T
```
`np.vstack()`：配列を縦に結合する
`.T`：行列の転置
#### yとXを使って、θを求める
```python
best_theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y.T)
# -> array([7.0962325 , 0.27487161])
```
`np.linalg.inv()`：逆行列
`A.dot(B)`：行列Aと行列Bの掛け算を行う
#### 等高線に描画
```python
plt.contour(theta_1_ax.T, theta_0_ax.T, Z, 100, cmap='jet')
plt.plot(theta_1_hist, theta_0_hist, 'x-')
# 正規方程式の結果を追加
plt.plot(best_theta[1], best_theta[0], 'x', color='red')
plt.xlabel('theta 1')
plt.ylabel('theta 0')
```
※赤が正規方程式で求めた$\theta$
![normal_and_gradient_descent_graph.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/8c535feb-1649-87c0-ebc0-b5bd61dfd882.jpeg)
**この結果から最急降下法で求めた$\theta$は最適解ではなかったが、徐々に近づいていたということがわかる。**-> あと少しだった...！
### ▶正規方程式の注意点
- 特徴量の数が多すぎると$(X^TX)^{-1}$の計算に時間がかかる
大体nが10000を超えるような高次のデータになってくると厳しい
⇒最急降下法を使用する
- 特徴量同士に強い相関がある(多重共線性：multicollinearity)と$(X^TX)^{-1}$を求めることができない
例：「駅からの距離」と「駅からかかった時間」
⇒相関がある特徴量はどちらかを落とす/組み合わせて新しい特徴量を作成する
- 特徴量数がデータ数を上回ると、回帰のパラメータを特定できない
⇒特徴量を減らす/データを増やす
### ▶ライブラリを使った線形回帰(Python)
この記事シリーズで扱う機械学習モデルは大体`scikit-learn`というライブラリに実装されている。
`import sklearn`で使用することができる。
- 線形回帰では`sklearn.linear_model.LinearRegression`クラスを使用する。
1．インスタンス生成
2．`.fit(X, y)`メソッドで学習(Xは行列、またはDataFrameの形である必要がある)
3．`.predicct(X)`で予測
4．`.coef_`で特徴量($\theta_1、\theta_2...$)、`.intercept_`で切片($\theta_0$)を取得可能

上記の1～4はどの機械学習モデルでも同じなので、流れを覚えておくとよい。
```python
from sklearn.linear_model import LinearRegression

# インスタンス生成
model = LinearRegression()

# データ準備
X = df['space'].values.reshape(-1, 1)
y = df['rent'].values

# 学習
model.fit(X, y)

# モデルの予測(入力は学習時と同じく二次元のArrayになることに注意)
model.predict([[30]]) # -> array([15.34238078])

# coefficient 係数
model.coef_ # -> array([0.27487161])

# intercept 切片
model.intercept_ # -> 7.096232499908833
```
## ■特徴量スケーリング(feature scaling)
### ▶特徴量スケーリングとは
- 前処理の一つで、尺度が異なる特徴量をそろえたりすることをいう
- 今回は特にメジャーな「標準化」と「正規化」を扱っていく
- 特徴量スケーリングはアルゴリズムによって必要だったり不要だったりする
#### 特徴量スケーリングが必要なアルゴリズム
最急降下法のような勾配法を使用したものや、距離を扱うもの

$\theta_j:=\theta_j-\alpha\frac{\partial}{\partial\theta_j}L(\theta)$
例えば上記のような最急降下法のパラメータ更新の式で、$\theta$の値が1000だったり、0.01だったりすると、$\theta$の更新スピードがそれぞれ異なり、不安定になる。


#### 特徴量スケーリングが不要なアルゴリズム
正規方程式のような解析的(数式を使って1発で解を出す)に解を求めるようなものや、決定木を使うもの
### ▶標準化(standardize)
- それぞれの特徴量空間のスケールに合わせる
- 平均0、標準偏差1にそろえる(平均を引いて標準偏差で割る)
※詳しくは[こちら](https://qiita.com/Yushin-Tati/items/af085499c49d244068ee#%E6%A8%99%E6%BA%96%E5%8C%96standardize)を参照

例：広さ(25~1000$m^2$)、築年数(0~50年)を特徴量とするとスケールの差が大きすぎる
### ▶Pythonで標準化
- `sklearn.preprocessiong.StandardScaler`クラスを使用する
１.インスタンス生成
２.`.fit(X)`で学習
３.`.transform(X)`で変換

:::note info
以下のコードで使用しているデータは下記gitに格納されています。
https://github.com/Yushin-Tati/Learning_machine_learning
・sample_data2.csv
:::
```python
# データのロード
import pandas as pd
df = pd.read_csv('sample_data2.csv')
X = df[['space', 'age']]

# 標準化
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()

scaler.fit(X)
X_scaled = scaler.transform(X)

# fit_transofrmで同時にfit&transformが可能
# X_scaled = scaler.fit_transform(X)
```
```text
・標準化後のデータ
array([[-1.13466215, -0.8062772 ],
       [-0.9410338 ,  0.63350351],
       [-0.82485679,  1.35339387],
       [-0.74740544,  0.2015693 ],
       [-0.51505142, -0.37434298],
       [ 0.06583364, -1.09423334],
       [ 0.22073632, -0.51832106],
       [ 0.49181602, -0.95025527],
       [ 1.3050551 , -0.51832106],
       [ 2.07956852,  2.07328422]])
```
```python
# Dataframe化
X_scaled_df = pd.DataFrame(X_scaled, columns=X.columns)

# 線形回帰
from sklearn.linear_model import LinearRegression
model = LinearRegression()
y = df['rent']
model.fit(X_scaled_df, y)
```
### ▶正規化(normalization)
- 値の範囲を0～1にrescaleする処理 ⇒ $\frac{x-x_{min}}{x_{max}-x_{min}}$
- 最小の値$x_{min}$は0に、最大の値$x_{max}$は1になる
よって、外れ値の影響を受けやすい
### ▶Pythonで正規化
- `sklearn.preprocessing.MinMaxScaler`クラスを使用
１.インスタンス生成
２.`.fit(X)`で学習
３.`.transform(X)`で変換
```python
from sklearn.preprocessing import MinMaxScaler
scaler = MinMaxScaler()
scaler.fit(X)
X_scaled = scaler.transform(X)
# 短縮型
# X_scaled = scaler.fit_transform(X)

# Dataframeに変換
pd.DataFrame(X_scaled, columns=X.columns)
```
||space|age|
|:-:|:-:|:-:|
|0|0.000000|0.090909|
|1|0.060241|0.545455|
|2|0.096386|0.772727|
|3|0.120482|0.409091|
|4|0.192771|0.227273|
|5|0.373494|0.000000|
|6|0.421687|0.181818|
|7|0.506024|0.045455|
|8|0.759036|0.181818|
|9|1.000000|1.000000|
### ▶特徴量スケーリングの注意点
- リスケールするのはあくまで特徴量毎で、データごとの特徴量ではない
- 目的変数(y)のリスケールは不要
- 学習データとテストデータを分ける場合は、学習データのリスケールに使ったパラメータ(**平均値や標準偏差**)をテストデータに使用する

![scaled_warning.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/a27e0f61-c7c5-2d1e-8f50-69d851b4a0a8.jpeg)
### ▶標準化vs正規化
- どちらのスケーリングを使用すればいいかの明確な基準はない
- 生のデータ、標準化したデータ、正規化したデータのそれぞれで結果を比較して決めるのが良い
- 標準化はデータが正規分布に従っている場合は特に有効、また外れ値の影響を受けたくない場合にも標準化を選ぶ
- 概ね、標準化のほうが選択されがち
## ■線形回帰の解釈
線形回帰は**真の関係($Y=f(X)+\epsilon$)が線形の形をしていると仮定して**行うが、実際に**線形に従うことの方が少ない**。そのため、他の非線形のモデルと比べると**精度は低いもの**になってしまう。

しかし、線形回帰はその**解釈性の高さ**からとても好まれている。
※形が限定されるので、**目的変数と特徴量の関係性を判断しやすい**。

実際の現場でも精度よりも解釈(関係性)が重要視されることの方が多い。
逆に言えば、精度が低い分、**関係性を導き出す必要がある**。
### ▶線形回帰の係数
- 特徴量の値が1増えた際の目的変数の増加量の平均値のこと
- 係数の大小はその特徴量の影響力の大小を表す
※特徴量毎の影響力なので、標準化していない限りは特徴量同士の**影響力を比較することはできない**
- 係数が0出ないことを統計的に判断することができる

![linear_regression_factor_example.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/c61fa283-e2f2-e482-bae5-a816d763b13d.jpeg)

:::note info
**なぜ係数が0でないことが重要なのか**
予測に関係ない特徴量(係数が0)のデータを無駄に取る必要がなくなり、モデルの軽量化/簡略化ができる。
:::
### ▶統計的仮設検定
必要であれば[こちら](https://qiita.com/Yushin-Tati/items/8aeb2c190e3494dd6b78)を参考に復習してください
### ▶線形回帰の係数の統計的仮設検定
- パラメータ$\theta$が0だと仮定して(帰無仮説)、t検定を行う
- 検定統計量tは自由度m-n-1のt分布に従う(m：データ数、n：特徴量数)
$t=\frac{\hat\theta_j-0}{SE(\hat\theta_j)}$ (SE:標準誤差、$\hat\theta_j$：j番目の特徴量のパラメータの推定値)
※標準誤差(standard Error)は推定量の標準偏差
▶ Pythonで係数の仮説検定(t検定)
- `statsmodels.api.OLS`クラスを使用
１.`OLS(y, X)`でインスタンス生成(ここのXはバイアス項$\theta_0$の欄も必要となるので注意)
２.`.fit()`で学習
３.`/summary()`で様々な統計定期情報を出力
- バイアス項の追加は`statsmodels.api.add_constant(X)`で可能
```python
# データ準備
import pandas as pd
df = pd.read_csv('sample_data2.csv')
y_col = 'rent'
X = df.drop(columns=[y_col])
y = df[y_col]

# 線形回帰
import statsmodels.api as sma
# バイアス項追加
X = sma.add_constant(X)
est = sma.OLS(y, X)
est_trained = est.fit()
print(est_trained.summary())
```
![sample_OLS.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/70c21561-17c6-eed6-ec38-6c4846d66de3.jpeg)
- coef：係数(coefficient)
- str err：SE($\theta$)
- t：t値
- P>|t|：p値(0.05以下なら有意差あり)

上記の表を見るとspaceはrentに関係しており、ageは関係していないように見える
ただ一般的に考えて、築年数が賃料に関係していないとは言い難いので、その場合はageがどのようなデータとなっているのかグラフにして確認する必要がある。
```python
import seaborn as sns
sns.pairplot(df)
```
![OLS_pairplot(df).jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/f5fadc47-2a06-0288-9aaf-b6e671c59610.jpeg)

:::note info
上図ではageのデータがばらけすぎていることがわかり、データが少ないか築年数の範囲が広すぎるのではという予想が立てられる
:::
**標準化した場合**
```python
# 標準化した場合
from sklearn.preprocessing import StandardScaler
X = df.drop(columns=[y_col])
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
X_scaled = sma.add_constant(X_scaled)
```
```text
X_scaled
array([[ 1.        , -1.13466215, -0.8062772 ],
       [ 1.        , -0.9410338 ,  0.63350351],
       [ 1.        , -0.82485679,  1.35339387],
       [ 1.        , -0.74740544,  0.2015693 ],
       [ 1.        , -0.51505142, -0.37434298],
       [ 1.        ,  0.06583364, -1.09423334],
       [ 1.        ,  0.22073632, -0.51832106],
       [ 1.        ,  0.49181602, -0.95025527],
       [ 1.        ,  1.3050551 , -0.51832106],
       [ 1.        ,  2.07956852,  2.07328422]])
```
```python
est = sma.OLS(y, X_scaled)
est_trained = est.fit()
print(est_trained.summary())
```
![OLS_pairplot(df)_standard.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/288a17e5-539b-b3a2-157a-ec829ae11365.jpeg)

:::note info
標準化してもt値やp値は変化していない
:::
### ▶線形回帰の係数の比較
- |$\hat\theta_j$|が大きい(影響力が大きい)と、|t|値は大きくなり、p値は小さくなる
→棄却域に入りやすくなる
- SE($\hat\theta_j$)が大きい(ばらつきが大きい)と、|t|値は小さくなり、p値は大きくなる
- 標準化した場合SE($\hat\theta_j$)がすべての$\hat\theta_j$で等しくなる
→それぞれの係数を比較することができる(重要度の比較)
### ▶係数のF検定
- 特徴量が一つでも関連があるかどうかを検定する
帰無仮説：$\theta_0=\theta_1=\theta_2=...=\theta_n=0$
対立仮説：少なくとも1つの$\theta$は0ではない
- F検定を行うことで第1種の誤りを防止することにつながる(t検定を行う前に実施する)
#### 復習：第1種の誤り、第2種の誤り
[こちら](https://qiita.com/Yushin-Tati/items/8aeb2c190e3494dd6b78#%E7%AC%AC1%E7%A8%AE%E3%81%AE%E8%AA%A4%E3%82%8A%E3%81%A8%E7%AC%AC2%E7%A8%AE%E3%81%AE%E8%AA%A4%E3%82%8A)参照
- 本当は目的変数に関連がない特徴量も「関連がある」と結論付けてしまう場合がある(第1種の誤り)
→特に特徴量が多い場合は注意

![Error_of_the_first_kind.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/bf846d5c-a744-edf0-c7eb-b6b5c5893387.jpeg)
→F検定を行うことでいずれかの特徴量に有意差があると出ても無視することができる
- F検定で使用する式は下記のようになる

![use_calculate_in_F_test.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/e7c689b5-00b6-c4dc-13b6-238b73e5fcce.jpeg)
### ▶検定統計量FとF分布
- 帰無仮説が正しい場合と対立仮説が正しい場合の残差の平方和を考える

例)特徴量が1つの場合
![RSS_example_F_test.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/9bb432ad-8f5c-2ef1-f400-89f3b8cd921a.jpeg)
- RSS$_0$とRSS$_a$の差の比率を考える
それぞれのRSSをそれぞれの自由度で割る

![RSSo_and_RSSa_F_test.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/81fe08fd-f07f-3717-00ee-10a3b7da04ed.jpeg)
#### 復習：F分布
[こちら](https://qiita.com/Yushin-Tati/items/8aeb2c190e3494dd6b78#f%E5%88%86%E5%B8%83)参照
- 検定統計量Fは自由度nおよび自由度m-n-1のF分布に従う
(誤差項$\epsilon$が正規分布に従うとき → そうじゃなくても近い分布にはなるので気にしなくてもいい)

![F_graph.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/195838ba-0771-32b2-445e-55370e3a27d5.jpeg)

:::note info
帰無仮説が正しいときF値は1となるので、F分布の1以上の範囲のp値を考えればいい
$\frac{TSS-RSS}{n}=\frac{RSS}{m-n-1}$
→$(m-n-1)(TSS-RSS)=nRSS$
→$(m-n-1)TSS=(m-1)RSS=1$ ($\because$n=0、TSS＝RSS)
:::
### ▶Pythonで係数の有意差(F検定)を確認
- t検定の時に出力した`.summary()`にすでにF値やp値について書かれている

![OLS_F_test.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/9a212aa0-2091-b850-1caa-c685c7ec7104.jpeg)
→ p<0.05なので、少なくとも1つは有意差があると分かる
#### t検定とF検定の関係
![F_test_result_relationship.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/e016015f-6906-1fe8-1338-a52296c6c7f0.jpeg)
## ■質的変数の特徴量
### ▶質的変数とは
- 「男女」や「曜日」など数値ではなく、カテゴリを値とする変数のこと
↔ 連続な値(量的変数、連続変数)
### ▶one-hotエンコーディング
- 質的変数を特徴量として使用する際のエンコーディング(線形回帰だけでなく、よく使われる)
- 質量変数をyes、no(1、0)の2次元行列に変換し、それぞれの特徴量と掛け合わせる
→ この際、質量変数の数がn個の場合、n×n-1の行列を作成する
- この時作成される変数は「ダミー変数」と呼ばれる

![explain_one_hot_encoding1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/d8447640-7a6a-7a72-7ac3-2ff53a96e4ea.jpeg)
![explain_one_hot_encoding2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/5ffd4b06-09b2-3ec3-e57c-3b2826bad5ad.jpeg)

:::note info
**なぜn-1列になるのか?**
行の値がすべて0の時、他の特徴量はすべて「否定」されている、ということになるので残った最後の特徴量が選択されているとわかるため。
:::
### ▶ダミー変数トラップ
one-hotエンコーディングをした際に、上図の「北」の列を作成してしまうとほかの特徴量と完全な相関関係を確立してしまう。⇒ **多重共線性(multicollinearity)**

2つの質的変数の特徴量を例にとると下記のようになる。($x_{i1}=1-x_{i2}$)
![dummy_trap_example.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/95c18919-2ea2-7303-9bb1-061017961097.jpeg)
### ▶Pythonでダミー変数作成
:::note info
下記git hubに格納されているテストデータを使用しています。
https://github.com/Yushin-Tati/Learning_machine_learning
sample_data3.csv
:::
- `pandas.get_dummies(df, drop_first=True)`で質的変数をone-hotエンコーディングできる
- `drop_first=True`で特徴量を1つ落としてくれる(ダミー変数トラップ対策)
- ダミー変数に対する標準化をしてもいいが、**しないのが一般的**
- ライブラリではそもそもダミー変数に対応されていることが多いので、`drop_first=True`をしなくても問題なく学習することができる
```python
import pandas as pd
df = pd.read_csv('sample_data3.csv')

# ダミー変数トラップを回避するためにdrop_first=Trueを指定
pd.get_dummies(df, drop_first=True)
```
![get_dummies.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/3a798acb-7bc6-a021-c4c4-82d7deafa2d4.jpeg)
#### 次の記事
[Python初心者の備忘録 #16 ~機械学習入門編02~](https://qiita.com/Yushin-Tati/items/5fe645f0a3c70b4dfa1a)