Weighted Least Squares in scikit-learn

# Reference

* https://medium.com/micin-developers/decipher-github-lr-sw-40e519a13c0a
* https://jp.mathworks.com/help/curvefit/least-squares-fitting.html#bq_5kr9-3
* https://www.mathworks.com/help/curvefit/least-squares-fitting.html#bq_5kr9-3

# Preparation

```
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
plt.rcParams['font.size']=15

def plt_legend_out(frameon=True):
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, frameon=frameon)
```

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

# OLS

ここでは、以下のデータを想定します。x=10の外れているデータ点を「外れ値」とみなすか、否かで対処が異なってきます。Ordinal Least Squaresモデルで無視するケース（WSL1）と組み込むケース（WSL2）を考えます。OLSですと、以下の通りx=10の点に引っ張られた回帰直線になります。

```
n = 11
x = np.arange(0,n,1)

np.random.seed(0)
y = x + np.random.randn(n)
y[len(y)-1] = y[len(y)-1]*2

plt.scatter(x,y,color='k')
plt.xlabel('x')
plt.ylabel('y')

reg = LinearRegression().fit(x.reshape(-1,1),y)
y_pr1 = reg.predict(x.reshape(-1,1))

plt.plot(x,y_pr1,label='pred')
plt_legend_out()
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/c270f8a3-7362-0fc4-ed88-86138abc07c5.png)

# WLS1

x=10の点を「外れ値」とみなして、処理を進めます。まずは、残差を取ります。

```
diff = y - y_pr1
df_plot = pd.DataFrame({'x':x,'error':diff})
sns.barplot(data=df_plot,x='x',y='error',color='k')
plt.ylabel('$y_i-\hat{y_i}$')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/c42b1703-e2b2-cbdd-7dec-5d00414f3d6e.png)

（略）

```
sw = (np.max(np.abs(diff)) - np.abs(diff))**2
df_plot = pd.DataFrame({'x':x,'error':sw})
sns.barplot(data=df_plot,x='x',y='error',color='k')
plt.ylabel('$(argmax|y_i-\hat{y_i}|-|y_i-\hat{y_i}|)^2$')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/baa9de0f-0daa-39f9-1c6f-369c42795991.png)

上記のsample weightをOrdinal Least Squaresに組み込みます。x=10の点を無視した回帰直線となりました。

```
reg.fit(x.reshape(-1,1),y, sample_weight=sw)
y_pr2 = reg.predict(x.reshape(-1,1))
plt.plot(x,y_pr1,label='Ordinal Least Squares')
plt.plot(x,y_pr2,label='Weighted Least Squares')
plt_legend_out()
plt.scatter(x,y,color='k')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/017fba11-91ef-a851-0fa7-a24797f54f29.png)

```
plt.figure(figsize=(10,4))
plt.subplots_adjust(wspace=0.4)

min = np.min(np.concatenate([y,y_pr1,y_pr2]))-1
max = np.max(np.concatenate([y,y_pr1,y_pr2]))+1

plt.subplot(1,2,1)
plt.scatter(y,y_pr1,color='k')
plt.plot([min,max],[min,max],color='gray',lw=0.5)
plt.xlabel('exp')
plt.ylabel('pred')
plt.grid()
plt.xlim(min,max)
plt.ylim(min,max)
plt.title('Ordinal Least Squares')

plt.subplot(1,2,2)
plt.scatter(y,y_pr2,color='k')
plt.plot([min,max],[min,max],color='gray',lw=0.5)
plt.xlabel('exp')
plt.ylabel('pred')
plt.grid()
plt.xlim(min,max)
plt.ylim(min,max)
plt.title('Weighted Least Squares')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/3fd80e65-b580-e1d8-edcb-3c9da24cce33.png)

# WLS2

x=10の点を過剰にフィットさせたいケースを考えます。

```
df_plot = pd.DataFrame({'x':x,'y':y})
sns.barplot(data=df_plot,x='x',y='y',color='k')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/14415fc8-4a5f-f230-832d-457a134ccd3c.png)

```
reg.fit(x.reshape(-1,1),y, sample_weight=1/y)
y_pr2 = reg.predict(x.reshape(-1,1))
plt.plot(x,y_pr1,label='Ordinal Least Squares')
plt.plot(x,y_pr2,label='Weighted Least Squares')
plt_legend_out()
plt.scatter(x,y,color='k')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/3fc6b70f-c69e-a836-9cfa-85c7d0ea97e1.png)

```
plt.figure(figsize=(10,4))
plt.subplots_adjust(wspace=0.4)

min = np.min(np.concatenate([y,y_pr1,y_pr2]))-1
max = np.max(np.concatenate([y,y_pr1,y_pr2]))+1

plt.subplot(1,2,1)
plt.scatter(y,y_pr1,color='k')
plt.plot([min,max],[min,max],color='gray',lw=0.5)
plt.xlabel('exp')
plt.ylabel('pred')
plt.grid()
plt.xlim(min,max)
plt.ylim(min,max)
plt.title('Ordinal Least Squares')

plt.subplot(1,2,2)
plt.scatter(y,y_pr2,color='k')
plt.plot([min,max],[min,max],color='gray',lw=0.5)
plt.xlabel('exp')
plt.ylabel('pred')
plt.grid()
plt.xlim(min,max)
plt.ylim(min,max)
plt.title('Weighted Least Squares')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/d3973b5e-4fe2-b66a-e4b8-a8cf02436da2.png)

さらにWeightを強く掛けます。

```
diff = y - y_pr1
df_plot = pd.DataFrame({'x':x,'error':np.abs(diff)*y})
sns.barplot(data=df_plot,x='x',y='error',color='k')
plt.ylabel('$y_i\cdot|y_i-\hat{y_i}|$')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/d5d3b787-4e54-32b2-5c01-54173da57974.png)

```
reg.fit(x.reshape(-1,1),y, sample_weight=np.abs(diff)*y)
y_pr2 = reg.predict(x.reshape(-1,1))
plt.plot(x,y_pr1,label='Ordinal Least Squares')
plt.plot(x,y_pr2,label='Weighted Least Squares')
plt_legend_out()
plt.scatter(x,y,color='k')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/031632ea-0115-1da0-9913-e3ec24508224.png)

```
plt.figure(figsize=(10,4))
plt.subplots_adjust(wspace=0.4)

min = np.min(np.concatenate([y,y_pr1,y_pr2]))-1
max = np.max(np.concatenate([y,y_pr1,y_pr2]))+1

plt.subplot(1,2,1)
plt.scatter(y,y_pr1,color='k')
plt.plot([min,max],[min,max],color='gray',lw=0.5)
plt.xlabel('exp')
plt.ylabel('pred')
plt.grid()
plt.xlim(min,max)
plt.ylim(min,max)
plt.title('Ordinal Least Squares')

plt.subplot(1,2,2)
plt.scatter(y,y_pr2,color='k')
plt.plot([min,max],[min,max],color='gray',lw=0.5)
plt.xlabel('exp')
plt.ylabel('pred')
plt.grid()
plt.xlim(min,max)
plt.ylim(min,max)
plt.title('Weighted Least Squares')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278468/3b199af9-1e53-3b88-9227-a652ccd0a00a.png)

# その他

```
>>> reg.predict(np.array(30).reshape(-1,1))
array([71.52016456])
```