教師あり学習の基礎その1ー単回帰ー（備忘録）

# 回帰
## （線形）単回帰モデル
1種類のデータから1種類の予測値を出力するモデルです。
入力データを$$x=(x_1, x_2, \cdots, x_n)$$、出力データを$$y=(y_1, y_2, \cdots, y_n)$$とするとき、__モデルを最もよく表している__直線$$y=wx+b$$の傾き$w$と切片$b$を決めます。
$w$を__重み__、$b$を__バイアス__といいます。
実際のデータとモデル（$y=wx+b$）の間には誤差があるので、その誤差を
$$\varepsilon=(\varepsilon_1, \varepsilon_2,\cdots , \varepsilon_n)$$とします。このとき、各教師データについて

```math
y_i = wx_i+b+\varepsilon_i
```

が成り立ちます。
モデルを最もよく表しているとは、適切な__損失関数__（例えば、二乗和誤差）を決め、モデル（直線の式）がそれを最小にするようなときにいいます。
損失関数$L(w,b)$を二乗和誤差（残差平方和ともいう）

```math
L(w, b):=\sum_{i=1}^n \{y_i-(wx_i+b)\}^2
```
にすると、$$L(w,b)=\sum_{i=1}^n\varepsilon_i^2$$すなわち誤差の二乗和となっていることが分かります。場合によっては、(1/2)倍します（微分すると係数が1にできるため）。
詳細は省きますが、$L(w,b)$を$w, b$それぞれで偏微分した式$=0$を解くことで、$w, b$を求めることができます。

```math
b = \bar{y}-w\bar{x}, \ w=\frac{Cov(x, y)}{Var(x)}=\frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x})}
```

ここで、$\bar{x}$は$x$の平均、$Var(x)$は$x$の分散、$Cov(x, y)$は$x$と$y$の共分散を表しています。
特に、この$(w, b)$の元で

```math
\begin{eqnarray}
&&y=wx+b=wx+(\bar{y}-w\bar{x})=\bar{y}+w(x-\bar{x})\\
&\Leftrightarrow& y-\bar{y}=w(x-\bar{x})
\end{eqnarray}
```

より、回帰直線は観測データの標本平均$(\bar{x}, \bar{y})$を通ります。

## 実装
Pythonで自力で実装＆sklearnのlinear_modelのLinearRegressionモデルを用いた実装をしてみます。
なお、再現性のために乱数シードを設定しています。
まずは必要なライブラリをインポート。

```python
import numpy as np
import matplotlib.pyplot as plt
from numpy.random  import randn
```

次に、線形回帰を関数として定義。

```python
def LinearRegression(x, y):
    n = len(x)
    temp_x = x.sum()/n
    temp_y = y.sum()/n
    w = ((x-temp_x)*(y-temp_y)).sum()/((x-temp_x)**2).sum()
    b = temp_y - w * temp_x
    return w, b
```

そして、教師データとして、$y=x+1+noise$を想定。$noise$は適当に標準正規分布に従う乱数に0.1をかけたものとしました。

```python
# initial value
# y=x+1にノイズを加えたものを想定
x = np.array([1,2,3,4,5])
np.random.seed(seed=0)
y = np.array([2,3,4,5,6])+np.random.randn(5)*0.1
print(y)
w, b = LinearRegression(x, y) 
print(w, b)
```

これを実行すると、およそ$w=1.02, b=1.08$と求まります。
念のため、グラフで可視化します。

```python
plt.scatter(x, y, color="k")
xx = np.linspace(0, 5, 50)
yy = xx * w + b
plt.plot(xx, yy)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/606074/35836cd5-3933-13fa-7225-d587b20b433b.png)

いい感じですね！

### ちょっとコードを使いやすいように書き換える
クラスを使ってコードを書き換えて見ます。
modelというクラスを作り、パラメータを求めるfitメソッド、予測を出力するpredictメソッド、決定係数を計算するscoreメソッドを作りました。決定係数$R^2$は、1に近いほど精度がよいということです。

```math
R^2:= 1-\frac{\sum_{i=1}^n (y_i-(wx_i+b))^2}{\sum_{i=1}^n (y_i-\bar{y})^2}
```

と定義されます。

```python
class model():
        
    def fit(self, X, y):
        n = len(X)
        temp_X = X.sum()/n
        temp_y = y.sum()/n
        self.w = ((X-temp_X)*(y-temp_y)).sum()/((X-temp_X)**2).sum()
        self.b = temp_y - self.w * temp_X
        
    def predict(self, X):
        return self.w*X +self.b
    
    def score(self, X, y):
        n = len(X)
        return  1 - ((y-(self.w*X+self.b))**2).sum() / ((y-y.sum()/n)**2).sum()
```

先ほどと同じように、教師データを入力します。
加えて、テストもさせてみます。

```python
X = np.array([1,2,3,4,5])
np.random.seed(seed=0)
y = np.array([2,3,4,5,6])+np.random.randn(5)*0.1

lr = model()
lr.fit(X, y)
print("w, b={}, {}".format(w, b))

# テストデータ
test_X = np.array([6,7,8,9,10])
pre = lr.predict(test_X)
print(pre)

# 決定係数R^2
print(lr.score(X, y))
```

グラフ化します。

```python
# 結果をグラフに表示
plt.scatter(X, y, color="k")
plt.scatter(test_X, pre, color="r")
xx = np.linspace(0, 10, 50)
yy = xx * w + b
plt.plot(xx, yy)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/606074/a912d9e5-7c58-59e8-ff66-27ecb75acb72.png)


## sklearnのライブラリを利用する
最後に、linear_modelの線形回帰LinearRegression()を用いて簡単にコードを書いてみましょう。

```python
from sklearn import linear_model

X = np.array([1,2,3,4,5])
np.random.seed(seed=0)
y = np.array([2,3,4,5,6])+np.random.randn(5)*0.1

X = X.reshape(-1,1)
model = linear_model.LinearRegression()
model.fit(X, y)
print(model.coef_[0])
print(model.intercept_)
# 決定係数R^2
print(model.score(X, y))
```

$w$はcoef_[0]、$b$はintercept_に相当します。値が一致していることがわかると思います。ここで注意する点は、1次元の入力データの場合は$X=np.array([[1],[2],[3],[4],[5]])$のように入力する必要があります。$X=np.array([1,2,3,4,5])$だとエラーを吐きます。
次に、先ほどのテストデータで予測をしてみます。

```python
test_X = test_X.reshape(-1,1)
pre = model.predict(test_X)
```

グラフ化します。

```python
plt.scatter(X, y)
plt.scatter(test_X, pre)
plt.plot(xx, model.coef_[0]*xx+model.intercept_)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/606074/ffb4c5f6-073c-cce9-2a92-37d2ee0192c9.png)

## まとめ
sklearnを用いてコードを書くと、簡単に次のようにかけます。これは、他の手法でもほぼ同様の形でかけます。あとはパラメータであったり、教師データを訓練データとテストデータに分けるなどオプション的な要素が加わるだけです。

```python
# モデルを定義
model = linear_model.LinearRegression()
# 学習(X:入力データ、y:正解ラベル)
model.fit(X, y)
# 予測(XX:未知のデータ)
model.predict(XX)
```

単回帰（最小二乗法）は以上です。
次回に続く。。。

## 追記
よく考えると、$b=(b, b, \cdots, b)\in \mathbf{\mathrm{R}}^n$とすべきでしたね。Pythonコードではブロードキャストのおかげでスカラー$b$としてもベクトルや行列との和は自動的に計算できます。

## 追記その2

```python
from sklearn import linear_model
```
としてますが、

```python
from sklearn.linear_model import LinearRegression
```

の方が楽ですね。他の書籍でもこちらが一般的かと。

## 参考文献

* 加藤公一『機械学習のエッセンス』SB Creative、2019年