Irisで前処理を試してみる。　～標準化、無相関化、白色化～


# はじめに
この記事では「[はじめてのパターン認識](https://www.morikita.co.jp/books/mid/084971)」(通称はじぱた)の第4章4.1節のデータの前処理(線形変換)の数値例をPython(numpyのみ)で実装する。


# データのロード
Scikit-kearnのload_iris()を使う。

```python
from sklearn.preprocessing import StandardScaler
import numpy as np
from matplotlib import pyplot as plt
```


```python
from sklearn.datasets import load_iris
from sklearn import preprocessing
iris = load_iris()
```


```python
import pandas as pd

df = pd.DataFrame(iris.data, columns=iris.feature_names)
df['target'] = iris.target
df.loc[df['target'] == 0, 'target'] = "setosa"
df.loc[df['target'] == 1, 'target'] = "versicolor"
df.loc[df['target'] == 2, 'target'] = "virginica"
```

データの中身を確認する。


```python
df.describe()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>sepal length (cm)</th>
      <th>sepal width (cm)</th>
      <th>petal length (cm)</th>
      <th>petal width (cm)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>count</th>
      <td>150.000000</td>
      <td>150.000000</td>
      <td>150.000000</td>
      <td>150.000000</td>
    </tr>
    <tr>
      <th>mean</th>
      <td>5.843333</td>
      <td>3.057333</td>
      <td>3.758000</td>
      <td>1.199333</td>
    </tr>
    <tr>
      <th>std</th>
      <td>0.828066</td>
      <td>0.435866</td>
      <td>1.765298</td>
      <td>0.762238</td>
    </tr>
    <tr>
      <th>min</th>
      <td>4.300000</td>
      <td>2.000000</td>
      <td>1.000000</td>
      <td>0.100000</td>
    </tr>
    <tr>
      <th>25%</th>
      <td>5.100000</td>
      <td>2.800000</td>
      <td>1.600000</td>
      <td>0.300000</td>
    </tr>
    <tr>
      <th>50%</th>
      <td>5.800000</td>
      <td>3.000000</td>
      <td>4.350000</td>
      <td>1.300000</td>
    </tr>
    <tr>
      <th>75%</th>
      <td>6.400000</td>
      <td>3.300000</td>
      <td>5.100000</td>
      <td>1.800000</td>
    </tr>
    <tr>
      <th>max</th>
      <td>7.900000</td>
      <td>4.400000</td>
      <td>6.900000</td>
      <td>2.500000</td>
    </tr>
  </tbody>
</table>
</div>



データの統計量を確認する。


```python
df.cov()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>sepal length (cm)</th>
      <th>sepal width (cm)</th>
      <th>petal length (cm)</th>
      <th>petal width (cm)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>sepal length (cm)</th>
      <td>0.685694</td>
      <td>-0.042434</td>
      <td>1.274315</td>
      <td>0.516271</td>
    </tr>
    <tr>
      <th>sepal width (cm)</th>
      <td>-0.042434</td>
      <td>0.189979</td>
      <td>-0.329656</td>
      <td>-0.121639</td>
    </tr>
    <tr>
      <th>petal length (cm)</th>
      <td>1.274315</td>
      <td>-0.329656</td>
      <td>3.116278</td>
      <td>1.295609</td>
    </tr>
    <tr>
      <th>petal width (cm)</th>
      <td>0.516271</td>
      <td>-0.121639</td>
      <td>1.295609</td>
      <td>0.581006</td>
    </tr>
  </tbody>
</table>
</div>



今回は「petal length」と「petal width」を対象にする。
```python
x=df['petal length (cm)']
y=df['petal width (cm)']
z=df['target']
```


```python
setosa = list(z.reset_index().query('target == "setosa"').index)
versicolor = list(z.reset_index().query('target == "versicolor"').index)
virginica = list(z.reset_index().query('target == "virginica"').index)
```

グラフ描画用の関数を作っておく。
```python
def data_plot():
    plt.rcParams['font.size'] = 14
    plt.rcParams['font.family'] = 'Times New Roman'

    #目盛を内側にする。
    plt.rcParams['xtick.direction'] = 'in'
    plt.rcParams['ytick.direction'] = 'in'

    #グラフの上下左右に目盛線を付ける。
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax1.yaxis.set_ticks_position('both')
    ax1.xaxis.set_ticks_position('both')

    #軸のラベルを設定する。
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    
    #元データ
    ax1.scatter(setosa_x, setosa_y, label='setosa', c='red')
    ax1.scatter(versicolor_x, versicolor_y, label='versicolor', c='blue')
    ax1.scatter(virginica_x, virginica_y, label='virginica', c='green')
    
    return ax1
```

```python
ax1=data_plot()

fig.tight_layout()
ax1.legend(loc="lower center", bbox_to_anchor=(0.5, 1.02,), borderaxespad=0, ncol=2)
```
![output_13_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/560543/28de951f-0200-6682-7875-fd9689f67a4e.png)


# 標準化
標準化とは測定単位の影響を取り除くべく、観測データの平均を0、分散を1にする変換のことである。
元のデータ$x$の平均$E\\{x\\}$を$\mu$、分散$Var\\{x\\}$を$\sigma$としたとき、線形変換
$$z=\frac{x-\mu}{\sigma}$$
を考えると、
$$E\\{z\\}=0、Var\\{z\\}=1$$
となる。
これらの関係が成り立つことを実装して確認する。

```python
# numpyで平均、分散それぞれを計算する。
# 平均
ave_x=np.average(x)
ave_y=np.average(y)
# 分散
stdx=np.std(x)
stdy=np.std(y)
```

    


```python
# 標準化の線形変換
nrom_x=[]
norm_y=[]
for i in range(len(x)):
    norm_x.append((x[i]-ave_x)/stdx)
    norm_y.append((y[i]-ave_y)/stdy)
```


```python
#setosa, viersicolor, virginicaそれぞれをプロットするために配列を準備
setosa_x=[]
versicolor_x=[]
virginica_x=[]
setosa_y=[]
versicolor_y=[]
virginica_y=[]
n_setosa_x=[]
n_versicolor_x=[]
n_virginica_x=[]
n_setosa_y=[]
n_versicolor_y=[]
n_virginica_y=[]

for i in setosa:
    setosa_x.append(x[i])
    setosa_y.append(y[i])
    n_setosa_x.append(norm_x[i])
    n_setosa_y.append(norm_y[i])

for i in versicolor:
    versicolor_x.append(x[i])
    versicolor_y.append(y[i])
    n_versicolor_x.append(norm_x[i])
    n_versicolor_y.append(norm_y[i])

for i in virginica:
    virginica_x.append(x[i])
    virginica_y.append(y[i])
    n_virginica_x.append(norm_x[i])
    n_virginica_y.append(norm_y[i])
```

【参考】Scikit-learnを使う場合


```python
# x,yを結合
data = np.c_[x, y]                                          

# 標準化を行う
sc = StandardScaler()
std = sc.fit_transform(data)
```

    

グラフを描画する。
```python
ax1=data_plot()

ax1.scatter(n_setosa_x, n_setosa_y, label='Standardized_setosa', marker='^', c='red')
ax1.scatter(n_versicolor_x, n_versicolor_y, label='Standardized_versicolor', marker='^', c='blue')
ax1.scatter(n_virginica_x, n_virginica_y, label='Standardized_virginica', marker='^', c='green')

fig.tight_layout()
ax1.legend(loc="lower center", bbox_to_anchor=(0.5, 1.02,), borderaxespad=0, ncol=2)
```

![output_21_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/560543/f4426c64-bb4c-5cda-a0f7-2b373d6f2a85.png)

  


 

# 無相関化
無相関化とは、観測データ間の相関をなくす処理のことである。観測データから得られた共分散行列$\Sigma$の固有値問題
$$\Sigma \boldsymbol{s}=\lambda \boldsymbol{s}$$
を解いて得られた固有値を$\lambda_{1} \geqq \lambda_{2} \geqq \cdot\cdot\cdot\geqq\lambda_{d} $、
対応する固有ベクトルを$\boldsymbol{s_{1}, s_{2}, \cdot\cdot\cdot, s_{d}}$として、行列$\boldsymbol{S}$を以下で定義する。

$$\boldsymbol{S}=(\boldsymbol{s_{1}, s_{2}, \cdot\cdot\cdot, s_{d}})$$

観測データ$\boldsymbol{x}$に対して$\boldsymbol{S}^{T}$で線形変換を行う。すなわち、$\boldsymbol{y}=\boldsymbol{S}^{T}\boldsymbol{x}$を考えたとき、以下が成り立つ。
$$E\\{\boldsymbol{y}\\}=\boldsymbol{S}^{T}\boldsymbol{\mu}, Var\\{\boldsymbol{y}\\}=\boldsymbol{S}^{-1} \boldsymbol{\Sigma}\boldsymbol{S}=\boldsymbol{\Lambda}$$

$\boldsymbol{S}^{-1} \boldsymbol{\Sigma}\boldsymbol{S}$は$\boldsymbol{\Sigma}$の対角化そのものであるため、対角成分以外は0となる。すなわち、共分散は0となり、相関係数の定義から相関係数は0となることがわかる。

```python
# 分散共分散行列の回転行列Sを求める。
_, S = np.linalg.eig(cov)           
S[:,1]=S[:,1]*(-1)
```
$\boldsymbol{S}^{T}=$$\boldsymbol{S}^{-1}$が成り立つ。すなわち、$S^{T}$$S$が単位行列であることが確認できる。
```python
np.dot(S.T, S)
```
    [[ 1.00000000e+00, -2.53792738e-17],
       [-2.53792738e-17,  1.00000000e+00]]

```python
# データを無相関化 ※.Tは転置を表す
pledged_decorr = np.dot(S.T, data.T).T
```

$\Lambda = \boldsymbol{S}^{-1} \boldsymbol{\Sigma}\boldsymbol{S}$を求める。
```python
lam = np.dot(np.dot(np.linalg.inv(S), cov), S)
```

はじぱたp.43の式(4.16)と同じ値であることを確認する。
```python
print(lam)
```
    [[3.66123805e+00 2.31134308e-16]
     [1.20962250e-16 3.60460707e-02]]


    


相関係数が0であることを確認する。
```python
print('相関係数: {:.3f}'.format(np.corrcoef(pledged_decorr[:, 0], pledged_decorr[:, 1])[0,1]))
```

    相関係数: 0.000
    

グラフを描画する。
```python
ax1=data_plot()

ax1.scatter(pledged_decorr[:, 0], pledged_decorr[:, 1], label='no correlation', marker='^', c='black')

plt.ylim(-3,3)

fig.tight_layout()
ax1.legend(loc="lower center", bbox_to_anchor=(0.5, 1.02,), borderaxespad=0, ncol=2)
```


![output_32_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/560543/d449c713-fcd7-d75c-4237-8669af0536ad.png)



# 白色化
前述の$\boldsymbol{\Lambda}、\boldsymbol{S}^{T}$、観測データ$\boldsymbol{x}$の平均$\boldsymbol{\mu}$を用いて、線形変換$\boldsymbol{u}=\boldsymbol{\Lambda}^{-\frac{1}{2}}\boldsymbol{S}^{T}(\boldsymbol{x}-\boldsymbol{u})$を考えると、以下が成り立つ。
$$E\\{\boldsymbol{u}\\}=0$$
$$Var\\{\boldsymbol{u}\\}=\boldsymbol{\Lambda}^{-\frac{1}{2}}\boldsymbol{S}^{-1}\boldsymbol{\Sigma}\boldsymbol{S}\boldsymbol{\Lambda}^{-\frac{T}{2}}
=\boldsymbol{\Lambda}^{-\frac{1}{2}}\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{\frac{1}{2}}=\boldsymbol{I}$$
ただし、$\boldsymbol{\Lambda}^{-\frac{T}{2}}$は$\boldsymbol{\Lambda}^{-\frac{1}{2}}$の転置行列である。このように観測データの標準偏差を1に正規化し、かつ中心化する変換を白色化と呼ぶ。

データを白色化する。
```python
ave=np.c_[ave_x, ave_y]
lam_1_2=np.linalg.inv(np.sqrt(lam))
u = np.dot(lam_1_2, np.dot(S.T, (df_corr-ave).T))
```

グラフを描画する。
```python
ax1=data_plot()
plt.xlim(-4,8)

ax1.scatter(u.T[:, 0], u.T[:, 1], label='Whitening', marker='^', c='black')

fig.tight_layout()
ax1.legend(loc="lower center", bbox_to_anchor=(0.5, 1.02,), borderaxespad=0, ncol=2)
```

![output_35_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/560543/7f184f5d-d5cb-2605-a3c5-78722d589fe7.png)


    


最後に白色化後のデータの共分散行列が単位行列になっていることを確認する。
```python
np.cov(u)
```
確かに単位行列になっている。



    array([[ 1.00000000e+00, -6.58744352e-08],
           [-6.58744352e-08,  1.00000000e+00]])