Pythonで多次元尺度法 (MDS) 〜距離行列から位置関係を再現する〜

n個体間の非類似度または距離が与えられているとき、それらn個体の位置関係を（低次元の）座標で表現する手法として、**多次元尺度法** (MDS : Multi-Dimensional Scaling) があります。

MDSの数理的な解説は別の機会に譲るとして、今回はscikit-learnのパッケージを使ってMDSを試してみます。MDSには大きく分けて計量MDSと非計量MDSに分けられますが、今回扱うのは計量MDSになります。

## ライブラリのインポート

scikit-learnでは `sklearn.manifold.MDS` をインポートすることでMDSサブパッケージを利用できます。

```python:In
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from tqdm import tqdm  # forループの進捗バーを表示するライブラリ
from sklearn.manifold import MDS
```

## データの読込み、前処理

[駅データ.jp](http://www.ekidata.jp/)より、`station20200619free.csv` をダウンロードし、pandas DataFrameとして読み込みます。

```python:In
data = pd.read_csv('station20200619free.csv')
print(data.head())
```

```:Out
   station_cd  station_g_cd station_name  station_name_k  station_name_r  line_cd  pref_cd      post           address         lon        lat    open_ymd   close_ymd  e_status   e_sort
0     1110101       1110101           函館             NaN             NaN    11101        1  040-0063    北海道函館市若松町１２-１３  140.726413  41.773709  1902-12-10  0000-00-00         0  1110101
1     1110102       1110102          五稜郭             NaN             NaN    11101        1  041-0813           函館市亀田本町  140.733539  41.803557  0000-00-00  0000-00-00         0  1110102
2     1110103       1110103           桔梗             NaN             NaN    11101        1  041-0801  北海道函館市桔梗３丁目４１-３６  140.722952  41.846457  1902-12-10  0000-00-00         0  1110103
3     1110104       1110104          大中山             NaN             NaN    11101        1  041-1121       亀田郡七飯町大字大中山  140.713580  41.864641  0000-00-00  0000-00-00         0  1110104
4     1110105       1110105           七飯             NaN             NaN    11101        1  041-1111         亀田郡七飯町字本町  140.688556  41.886971  0000-00-00  0000-00-00         0  1110105
```

緯度(lat)・経度(lon)の情報をもとにmatplotlibで描画してみます。

```python:In
fig, ax = plt.subplots()
data.plot('lon', 'lat', kind='scatter', ax=ax, marker='.')
plt.show()
```

![スクリーンショット 2020-12-04 15.58.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/4c9a8172-75d2-f5ee-eb46-966521e15e06.png)

日本ですね！ただこのままだとデータ数が多いので、九州部分だけを切り出してみます。

```python:In
data = data[(data['lon']>128)&(data['lon']<132)]
data = data[data['lat']<34]
len(data)
# 1207
```

```python:In
fig, ax = plt.subplots()
data.plot('lon', 'lat', kind='scatter', ax=ax, marker='.')
plt.show()
```

![スクリーンショット 2020-12-04 16.03.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/c81b699f-bc03-5b17-c23e-d1634b1045ce.png)


## 距離行列の計算

抽出したデータをもとに各駅間の距離行列$\boldsymbol{D}$を計算します。ユークリッド距離行列の場合、駅間距離の二乗の値になります。

```python:In
N = len(data)
D = np.zeros((N, N), dtype=float)
for i in tqdm(range(N)):
    for j in range(N):
        lon_diff = (data['lon'].iloc[i] - data['lon'].iloc[j])
        lat_diff = (data['lat'].iloc[i] - data['lat'].iloc[j])
        D[i][j] = lon_diff**2 + lat_diff**2

print(D)
```

```:Out
array([[0.00000000e+00, 1.26355528e-03, 6.53891820e-03, ..., 1.06203018e+00, 1.06368290e+00, 1.05709741e+00],
       [1.26355528e-03, 0.00000000e+00, 2.39660029e-03, ..., 9.92989742e-01, 9.94425737e-01, 9.87970532e-01],
       [6.53891820e-03, 2.39660029e-03, 0.00000000e+00, ..., 9.02220608e-01, 9.03831353e-01, 8.97815038e-01],
       ...,
       [1.06203018e+00, 9.92989742e-01, 9.02220608e-01, ..., 0.00000000e+00, 6.80062050e-05, 1.68960250e-04],
       [1.06368290e+00, 9.94425737e-01, 9.03831353e-01, ..., 6.80062050e-05, 0.00000000e+00, 3.12263650e-05],
       [1.05709741e+00, 9.87970532e-01, 8.97815038e-01, ..., 1.68960250e-04, 3.12263650e-05, 0.00000000e+00]])
```

## 低次元空間へのマッピング

距離行列$\boldsymbol{D}$のみを用いて、もとの二次元空間の位置関係を再現していきます。`sklearn.manifold.MDS` クラスに次の引数を与えてモデルを生成し、`fit_transform()` メソッドで適合・変換します。

* `n_components` ... マッピング先の次元数。今回は二次元空間なので2。
* `dissimilarity` ...
    - `euclidean` または `precomputed` (デフォルト : `euclidean`)
    - 距離行列を既に計算している場合は `precomputed` を指定。
* `random_state` ... 再現性を持たせるために乱数シードを設定。


```python:In
mds = MDS(n_components=2, dissimilarity="precomputed", random_state=0)
pos = mds.fit_transform(D)
```

## 結果の確認

```python:In
res = pd.DataFrame(pos, columns=['x', 'y'])
plt.scatter(res['x'], res['y'], marker='.')
```

![スクリーンショット 2020-12-04 15.59.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/ca8a57b7-83ae-56d5-729c-297ece6c627c.png)

九州のようなものが見えます。軸が回転しているので戻してみしょう。

```python:In
def inverse(x):
    return x*(-1)

plt.scatter(res2['y'].map(inverse), res2['x'], marker='.')
```

![スクリーンショット 2020-12-04 15.59.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/ad95d2a9-f370-c450-75e1-131f8a4beb5a.png)

ちょっと歪ですが九州ですね！ということで無事、距離行列のみから平面上の位置関係を再現することが出来ました。

## 参考

* [日本統計学会公式認定 統計検定準1級対応 統計学実践ワークブック](https://www.amazon.co.jp/dp/478060852X/)
* [[Python][Scikit-learn] 多次元尺度構成法(MDS)で都内公園の時間地図を作成する - Developers.IO](https://dev.classmethod.jp/articles/mds-for-parks/)
* [2.2.8. Multi-dimensional Scaling (MDS) - scikit-learn](https://scikit-learn.org/stable/modules/manifold.html#multidimensional-scaling)