k-NN法をフルスクラッチ実装

今日はk-NN法(k-NearestNeighbor)を実装する.

k-NN法の概念

k-NN法は教師データから未知データを予測する教師あり学習の一つで, また, 線形回帰などのようにパラメータを最適化するような手法を取らない.
すなわち, 「教師データから識別境界のようなものを学習してから未知データが境界のどちらにあるか」という手法をとらずに, 「それぞれの未知データに対して教師データとの関係を見てどちらのクラスかを予測する」という方法をとる. これを怠惰学習というらしい.
k-NN法に関しては, 名前の通り, 未知データに対して最も近い教師データk個のうち多い方のクラスラベルを予測クラスとする.
まず下図のようなデータがあるとする. 赤色がクラス1, 緑色がクラス0 とする. 星がクラスラベルのわかっていない未知データである.
この未知データがどちらのクラスに入れた方が良いのかを予測したい.

f:id:umashika5555:20181223013754p:plain
データ


今k=3すると未知データから最も近い3つの教師データを抽出できる.
f:id:umashika5555:20181223020457p:plain
3-NN
赤色(クラス1)のデータが1つで緑色(クラス0)のデータが2つである.
よって多数決をとって, 緑(クラス0)データの方が多いので, 未知データは緑(クラス0)と予想する.
式で表すとこのようになるが, わざわざ式で表さなくてもよい.
f:id:umashika5555:20181223015047g:plain

データの生成

k-NN法を実装するために二次元のデータを生成する.
k-NN法は分布を仮定した方法ではないため, 分布から生成する必要は無いと思うが, ラベル付きのデータを生成したかったので取り敢えずガウス分布を使ってデータを生成した.

f:id:umashika5555:20181223015313p:plain
データ

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 乱数固定
np.random.seed(0)
# 平均ベクトル
Mu1, Mu2 = np.array([1, 10]), np.array([9,2])
# 分散共分散行列
sigma11, sigma12, sigma21, sigma22 = 4, 5,  6, 3
SIGMA1, SIGMA2 = np.array([[sigma11**2, np.sqrt(sigma11*sigma12)],[np.sqrt(sigma11*sigma12), sigma12**2]]), np.array([[sigma21**2, np.sqrt(sigma21*sigma22)],[np.sqrt(sigma21*sigma22), sigma22**2]])
# データ生成
values1 = np.random.multivariate_normal(Mu1, SIGMA1, 200)
values2 = np.random.multivariate_normal(Mu2, SIGMA2, 200)
# データ保存
np.save("./data/values1.npy", values1)
np.save("./data/values2.npy", values2)
# 散乱図の表示
plt.scatter(values1[:,0], values1[:,1], color="r")
plt.scatter(values2[:,0], values2[:,1], color="g")
plt.xlabel("x")
plt.ylabel("y")
plt.xlim(-10, 30)
plt.ylim(-10, 25)
plt.savefig("./img/data.png")
plt.show()

k-NN法による未知データの分類

格子点を未知データとして考え領域分割を行う.
格子点から各データ点までの距離を計算し, 最も近いデータ点k個のラベルから多数決を行う.
閾値が0.5のとき, kが奇数の場合は多数決が成立するのだが, kが偶数の場合で各ラベルが同じ個数あるとき多数決が成立しない.
よってkは奇数にするほうが良いことが分かる.
またkが大きくなると, 汎用性が高くなることが分かる.

f:id:umashika5555:20181223015636p:plain
k=1
f:id:umashika5555:20181223015651p:plain
k=2
f:id:umashika5555:20181223015705p:plain
k=3
f:id:umashika5555:20181223015720p:plain
k=4
f:id:umashika5555:20181223015818p:plain
k=5

f:id:umashika5555:20181223015836p:plain
k=7

f:id:umashika5555:20181223015854p:plain
k=9

f:id:umashika5555:20181223151116p:plain
k=51

import numpy as np
import matplotlib.pyplot as plt

# k-NN法のパラメータkを設定
k = 200
# データのロード
values1 = np.load("./data/values1.npy")
values2 = np.load("./data/values2.npy")

# データを
# 統合, その際に教師情報もつける
data = np.vstack((values1, values2))
labels = np.hstack((np.ones(len(values1)), np.zeros(len(values2))))

# 未知のデータとして格子点を作成
# 2つのデータから最小, 最大のx, yを探す
x_min, x_max = min(np.min(values1[:,0]), np.min(values2[:,0])), max(np.max(values1[:,0]), np.max(values2[:,0]))
y_min, y_max = min(np.min(values1[:,1]), np.min(values2[:,1])), max(np.max(values1[:,1]), np.max(values2[:,1]))
# 格子点を作成
xx = np.linspace(x_min-1, x_max+1, 300)
yy = np.linspace(y_min-1, y_max+1, 300)
xxx, yyy = np.meshgrid(xx, yy)# 格子点(xxx, yyy)が生成された. これを未知データとみなしてk-NN法を適用する.

class0 = []
class1 = []
thresholds = []

# k-NN法の適用
for i, (xx, yy) in enumerate(zip(xxx, yyy)):
    for j, (x, y) in enumerate(zip(xx, yy)):
        # 未知データ(x,y)から既知データdata各点への距離を計算する
        # 距離はユークリッド距離の二乗, すなわちl2ノルムの二乗を計算する
        tmp = data - (x, y)
        distances = np.linalg.norm(tmp, axis=1)
        # 最も小さいスコア上位k点のラベルを参照
        supervisers_index = np.argsort(distances)[:k]
        # print(supervisers_index)
        supervisers = labels[supervisers_index]
        # print(supervisers)
        res = np.sum(supervisers)/ k 
        # if res > 0.5 then x_res in class1
        if res > 0.5:
            class1.append((x, y))
        elif res == 0.5:
            thresholds.append((x, y))
        else:
            class0.append((x, y))
        
class0, class1, thresholds = np.array(class0), np.array(class1), np.array(thresholds)
plt.scatter(values1[:,0], values1[:,1], color="r", marker="x", s=20)
plt.scatter(values2[:,0], values2[:,1], color="g", marker="x", s=20)
plt.scatter(class1[:,0], class1[:,1], color="r", marker="o", s=1, alpha=0.1)
plt.scatter(class0[:,0], class0[:,1], color="g", marker="o", s=1, alpha=0.1)
try:# 奇数の場合, 必ず多数決が成立するので閾値(0.5)とイコールになるものがないためエラー処理しておく
    plt.scatter(thresholds[:,0], thresholds[:,1], color="b", s=1, alpha=0.1)
except:
    pass

plt.title("./img/{}-NN methods".format(k))
plt.xlabel("x")
plt.ylabel("y")
plt.xlim(-9, 23)
plt.ylim(-7, 23)
plt.savefig("./img/{}-NN.png".format(k))
plt.show()

最小二乗法を実装する

予測モデルをつくるときに最も簡単な方法としては「最小二乗法による線形モデル」と「k近傍モデル」である.
この2つは統計学機械学習を学んでいない人でも割と思いつきそうな方法である.
線形モデルとはパラメータに対して線形なモデルという意味である.
例えば, y = ax + b という1次関数を中学で学んだ.
今, データxが観測されたときに, yを予測するというタスクを行う.
このときa(傾き), b(切片)の値によって予想するべき値yが異なってくる.
すなわちa, bが決定されれば全てのxの入力に対してyという出力ができるようになる.
このようなa, b をパラメータという.
これを多次元に拡張すると

y = w_0 + w_1x_1 + w_2x_2 + \ldots w_nx_n
これは

y = w_01 + w_1x_1 + w_2x_2 + \ldots w_px_p
と見て

y = \boldsymbol{x}^T \boldsymbol{w}
とベクトル表記で書くことができる.
これは, p+1 次元の一つのデータに対する予測の式なので, n個のデータを同時に扱う場合には

\boldsymbol{y} = \boldsymbol{X}^T \boldsymbol{w}
とする.  \boldsymbol{X} は(p+1, n)行列であり転置すると(n, p+1)行列, パラメータベクトル \boldsymbol{w}は(p+1)ベクトルである. ((p+1, 1)行列という表現もできる)
この線形モデルを, 最小二乗法によりフィットする. (予測と実際の誤差を最小にする)
 
min: RSS(\boldsymbol{w}) = (\boldsymbol{y} - \boldsymbol{X}\boldsymbol{w})^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{w})

\begin{align*}
RSS(\boldsymbol{w}) &= (\boldsymbol{y}^T-(\boldsymbol{X}\boldsymbol{w})^T)(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{w})\\
 &= \boldsymbol{y}^T\boldsymbol{y}-\boldsymbol{y}^T\boldsymbol{X}\boldsymbol{w}-(\boldsymbol{X}\boldsymbol{w})^T\boldsymbol{y}+(\boldsymbol{X}\boldsymbol{w})^T(\boldsymbol{X}\boldsymbol{w})\\
 &= \boldsymbol{y}^T\boldsymbol{y}-(\boldsymbol{y}^T\boldsymbol{X})\boldsymbol{w}-\boldsymbol{w}^T(\boldsymbol{X}^T\boldsymbol{y})+\boldsymbol{w}^T(\boldsymbol{X}^T\boldsymbol{X})\boldsymbol{w}
\end{align*}
この式はパラメータ\boldsymbol{w}に関して二次式であるので最小値が存在する.
パラメータ\boldsymbol{w}微分して\boldsymbol{0}になるパラメータ\boldsymbol{w}を求める.

\begin{align*}
&\boldsymbol{0} - \frac{\partial (\boldsymbol{y}^T\boldsymbol{X})\boldsymbol{w}}{\partial \boldsymbol{w}} - \frac{\boldsymbol{w}^T(\boldsymbol{X}^T\boldsymbol{y})}{\partial \boldsymbol{w}} + \frac{\partial \boldsymbol{w}^T(\boldsymbol{X}^T\boldsymbol{X})\boldsymbol{w}}{\partial \boldsymbol{w}} = \boldsymbol{0}\\
&-(\boldsymbol{y}^T\boldsymbol{X})^T-(\boldsymbol{X}^T\boldsymbol{y}) + (\boldsymbol{X}^T\boldsymbol{X} + (\boldsymbol{X}^T\boldsymbol{X})^T)\boldsymbol{w} = \boldsymbol{0}\\
&-\boldsymbol{X}^T\boldsymbol{w}-\boldsymbol{X}^T\boldsymbol{y} + (\boldsymbol{X}^T\boldsymbol{X} + \boldsymbol{X}^T\boldsymbol{X})\boldsymbol{w} = \boldsymbol{0}\\
&-2\boldsymbol{X}^T\boldsymbol{y} + 2\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{w} = \boldsymbol{0}\\
&2\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{w} = 2\boldsymbol{X}^T\boldsymbol{y}\\
&\boldsymbol{w} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}
\end{align*}


ベクトルに対する微分の以下の式は覚えておいたほうが良い.

\begin{align*}
\frac{\partial A^T \boldsymbol{x}}{\partial \boldsymbol{x}} = \frac{\partial \boldsymbol{x}^T A}{\partial \boldsymbol{x}} &= A\\
\frac{\partial \boldsymbol{x}^TA\boldsymbol{x}}{\partial \boldsymbol{x}} &= (A+A^T)\boldsymbol{x}\\
if\  A\  is\  symmetric\  matrix\\
&= 2A\boldsymbol{x}
\end{align*}

線形モデルの最小二乗法について, パラメータの解が

\begin{align*}
\boldsymbol{w} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}
\end{align*}
で得られることが分かった.
これをPythonによって実装してみる. 今回は一次関数のパラメータa,b を求め直線フィッティングを行う.
まずデータを生成する. データ生成のアルゴリズムは以下の通りに行った.

  1. a, b, σを決定する
  2. [-10, 10]の範囲でデータxをランダムに生成する
  3. 決定したa, bからt=ax+bを計算する
  4. t=ax+b+εを計算する. ただしε~N(0, σ^2)
### 単純な線形モデルy = ax + b についてのパラメータa,b を求める問題におけるデータを生成する
### 
### a, b, σの値を決定する
### [-x_min, x_max]のデータxをランダムにN個用意 
### 各データに対して, N(ax+b, σ^2)に従うガウス分布からtを生成
### 

import numpy as np
from numpy.random import seed
import matplotlib.pyplot as plt

seed(0)

a, b, sigma = 2, 3, 3
N = 100
x_min, x_max = -10, 10

# データ生成
X = (x_max-x_min) * np.random.rand(N) + x_min
Epsiron = np.random.normal(0, sigma**2, N)# N(0, sigma^2)のN個の正規乱数
T = (a*X + b) + Epsiron

# データをファイルに出力
f = open("data00.csv", "w")
for (x, t) in zip(X, T):
    f.write("{}, {}\n".format(x,t))
f.close()

# データをプロット
plt.scatter(X, T)
plt.show()

以下の図のようになった.

f:id:umashika5555:20181221091859p:plain
データの様子
このデータから,

\begin{align*}
\boldsymbol{w} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}
\end{align*}
に従ってパラメータを決定する.

import numpy as np
import matplotlib.pyplot as plt

# データの読み込み
name_data = "data00"
path_data = name_data + ".csv"
data = np.loadtxt(path_data,delimiter=",")
x, t = data[:,0], data[:,1]
p = x.shape[0]
ones = np.ones(p)
x_ = np.c_[ones, x]

# パラメータベクトルの計算
A = np.dot(x_.T, x_)
w = np.dot(np.dot(np.linalg.inv(A), x_.T), t)
print(w)

# データのプロット
plt.scatter(x, t)

# 回帰曲線の描画
b, a = w
xx = np.linspace(np.min(x)-2, np.max(x)+2, 100)
yy = a * xx + b
plt.plot(xx, yy, color="red")

plt.xlabel("x")
plt.ylabel("y")
plt.title("a=2, b=3, σ=2")
plt.savefig(name_data.format(a,b))
plt.show()

f:id:umashika5555:20181221092149p:plain
data00
上図のように実際に目で観ても綺麗にフィットできていることがわかる.


次回はこの線形モデルを更に発展させ, 基底関数, 正則化あたりをやっていきたい.

convert コマンド よく使う操作

画像のリサイズ

a.pngを指定した高さ/幅に同一比のまま縮尺拡大したb.pngを出力する.

# 幅に合わせる
convert -resize 640x a.png b.png
# 高さに合わせる
convert -resize x640 a.png b.png

a.pngという画像を640px * 640pxのb.pngに変換する.

convert -resize 640x640! a.png b.png

ディレクトリ内にある全ての画像に適応するならば

for filename in *.JPG; do convert -resize 640x640! ${filename} ${filename%.JPG}_640x640.jpg; done;

# ファイル名を変えずに実行したい場合
mogrify -resize 640x640! *.png;

ディレクトリresize_dirに変換後の画像を格納する場合

mkdir resize_dir
for filename in *.JPG; do convert -resize 640x640 $filename reseize_dir/${filename%.JPG}_640x640.JPG; done;

画像の形式変換

ディレクトリ内のJPGファイルをPNGファイルに一括変換する.

# ディレクトリ内のJPGファイルを探す
find ./ -name "*.jpg"
# ディレクトリ内のJPGファイルの個数を数える
find ./ -name "*.jpg" | wc -l
# JPGファイルをPNGに変換する
for filename in *.jpg;do convert "$filename" "${filename%.jpg}.png";done

【環境設定】ubuntu 18.04 LTSにslackを入れる

slackのデスクトップアプリケーションを, dpkgでインストールする方法
単純にslackのwebサイトからダウンロードして

$ sudo dpkg -i xxx.deb

とやったらlibappindicator1が無いみたいなエラーがなされたので,

$ sudo apt --fix-broken install
$ sudo apt install libindicator1
$ sudo dpkg -i xxx.deb

と行う.

【環境設定】ubuntu18.04のproxy関連設定

aptの設定

自宅でインストールしたubuntuをプロキシ環境下にあるネットワーク環境に持っていったとき, システム全体のプロキシ設定でブラウジングなどできたので, apt installなどできると思っていたが別に設定しないといけないらしい.
参考にしたのはこれ
usado.jp

$ sudo emacs /etc/environment

http_proxy="http://proxy-server:port/"
https_proxy="http://proxy-server:port/"

gitの設定

qiita.com
qiita.com

$git config --global http.proxy http://proxy.example.com:8080
$git config --global https.proxy http://proxy.example.com:8080

設定場所は.gitconfig

BitBucketの設定

22ポートが閉じられているので443ポートで通信を行う.

$ touch ~/.ssh/config
$ emacs ~/.ssh/config
Host bitbucket.org
    User git
    IdentityFile ~/.ssh/id_rsa
    HostName altssh.bitbucket.org
    Port 443
    ProxyCommand connect.exe -H hoge.proxy.jp:1080 %h %p

qiita.com
qiita.com
qiita.com
re-dzine.net

【環境設定】Ubuntu 18.04LTSでCtrl + Spaceで変換する

mozcを起動してデフォルトでHankaku/Zenkakuの項目をメモして
f:id:umashika5555:20180917024215p:plain
全く同じものをCtrl+Spaceに置き換えてエントリーを追加すればよい.
f:id:umashika5555:20180917024220p:plain
f:id:umashika5555:20180917024217p:plain
mozcを閉じるときに「設定は新しいアプリケーションから適応されます」的なことを言われたので, 新しいアプリを起動して試してみたができなかった.
しかし, 再起動したらできた.

【環境設定】ubuntuにGPUドライバを入れる

自分のノートPCにubuntu 18.04 LTSを入れて, GPUドライバがデフォルトでIntelの内蔵のグラフィックドライバになっていたため, NVIDIAのに変えた.
ここを参考にしたらスンナリできた.
taktak.jp

$ ubuntu-drivers devices
$ sudo ubuntu-drivers autoinstall
$ sudo reboot
$ nvidia-smi

【urllib】プロキシ設定

proxyの使い方メモ
keras.datasets.mnist.load_data()とかしたいときに, プロキシを通す方法

import urllib
proxy_support = urllib.request.ProxyHandler({'https': 'http://proxy.hogehoge.ac.jp:80'})
opener = urllib.request.build_opener(proxy_support)
urllib.request.install_opener(opener)

【Python】matplotlibで最低限のグラフをかけるようにする

TeXで表示させるため, グラフが必要な時期になってきた.
Pythonで最低限のグラフを描画する方法をメモ.
sigmoid関数のプロット

import numpy as np
import matplotlib.pyplot as plt

# 関数定義
x = np.linspace(-8, 8, 100)
y = 1/ (1+np.exp(-x))

# 垂直線
plt.vlines(0.0, -0.1, 1.1, colors="r", linestyle="dotted", label="")
# 水平線
plt.hlines(0.0, -8.0, 8.0, colors="k", linestyle="dotted", label="")
plt.hlines(0.5, -8.0, 8.0, colors="k", linestyle="dotted", label="")
plt.hlines(1.0, -8.0, 8.0, colors="k", linestyle="dotted", label="")

# sigmoid 関数
plt.plot(x, y, label="sigmoid function")

# y軸の設定
yticks = [0,0.5,1.0]
plt.yticks(yticks)

# 軸の名称
plt.xlabel(r"z")
plt.ylabel(r"$\phi(z)$")

# 凡例
plt.legend(loc="upper left")

# 図の保存
plt.savefig("./sigmoid.png")
plt.show()

f:id:umashika5555:20180122130503j:plain
TeXファイルに貼り付ける場合は, .epsにしておく必要がある場合がある.
convert コマンドでファイルの種類を変更する.

$ convert sigmoid.png sigmoid.eps

【Python】 pythonからshellの実行

Pythonからshellコマンドを実行する際に以下のようなやり方がある.
今回は次のshellコマンドをpythonで実行するやり方を比較する.

$sha256sum hoge.zip
#返り値:  xxxxxxxxxxx hoge.zip

os.system()を使う方法

import os
file_name = "hoge.zip"
os.system('sha256sum "{}"'.format(file_name))
# 成功
# xxxxxxxxxxx hoge.zip
# 0
# 失敗
# sha256sum: hoge.zip : No such file or directory
# 256

成功したか失敗したかのみしかわからない.
コマンドの出力が受け取れないし, 推奨されていないやり方らしい.

subprocess.check_output()を使う方法

import subprocess
file_name = "hoge.zip"
args = ["sha256sum", file_name]
subprocess.check_output(args)
# 成功
# b'16e3a87f14dce201f90dc0f7dd2a6318e7b9e5aef095e5815123c3a44e92ac97  GA_result.txt\n'
# 失敗
# sha256sum: hoge.zip : No such file or directory
#Traceback (most recent call last):
# File "<stdin>", line 1, in <module>
#  File ".pyenv/versions/3.5.4/lib/python3.5/subprocess.py", line 316, in check_output  **kwargs).stdout
#  File ".pyenv/versions/3.5.4/lib/python3.5/subprocess.py", line 398, in run output=stdout, stderr=stderr)
# subprocess.CalledProcessError: Command '['sha256sum', 'GA_result1.txt']' returned non-zero exit status 1

成功した場合, 返り値はbyte型なのでstr型にdecodeし, pythonコード内でshellの出力結果を得ることができる.

res = subprocess.check_output(args)
hash_value = res.decode("utf-8").split(" ")[0]

ただし, windowsの場合はうまく行かない可能性がある. 原因はまだ良く見てない.

commands.getstatusoutput()を使う方法

昔よく使っていたが, Python3 で廃止されたらしい.
これの代用がsubprocessモジュールになったらしい.

【TeX】 ページを跨いで表を表示する

longtable.styを使う

\newcolumntype{A}{>{\raggedright}p{0.3cm}}
\newcolumntype{N}{>{\raggedright}p{3.6cm}}
\newcolumntype{D}{>{\raggedright}p{10.0cm}}
{\footnotesize
\begin{longtable}[c]{|A|N|D|}
\hline
No & permission-name & Description \tabularnewline \hline
\endfirsthead
0  & \verb|ACCESS_CHECKIN_PROPERTIES| & Allows read/write access to the "properties" table in the checkin database, to change values that get uploaded.\tabularnewline \hline
1  & \verb|ACCESS_COARSE_LOCATION| & Allows an app to access approximate location. \tabularnewline \hline
2  & \verb|ACCESS_FINE_LOCATION| & Allows an app to access precise location.\tabularnewline \hline
3  & \verb|ACCESS_LOCATION_|\\ \verb|EXTRA_COMMANDS| & Allows an application to access extra location provider commands. \tabularnewline \hline
\end{longtable}
}

f:id:umashika5555:20180118020916p:plain

【参考】
http://www.biwako.shiga-u.ac.jp/sensei/kumazawa/tex/longtable.html
https://ctan.org/pkg/longtable

【Python】 再帰的にファイルを取得する

import os
def find_all_files(directory):
    for root, dirs, files in os.walk(directory):
        #yield root # ここをアンコメントするとディレクトリも得られる
        for file in files:
            yield os.path.join(root, file)
        
if __name__ == "__main__":
    path_directory = "SAMPLE"
    for file in find_all_files(path_directory):
        print(file)

treeコマンドのようにファイルを取得できる.
os.walk()とos.path.join()は始めて知った.
os.path.join(path_directory, name_file)でpath_directory + "/" + name_fileとしてくれるようだ.
普段

path ="aaaa"
for file in os.listdir(path):
    path_file = path + "/" + file

とファイルのpathを取得していたためこれは少し感動.
さらにwindowslinuxの"\\"と"/"の違いも一々気にする必要なさそうで便利!

参考
https://qiita.com/suin/items/cdef17e447ceeff6e79d

【Cygwin】 ホームディレクトリを変更する

Cygwinのterminal上ではなく, ホストOSであるWindows環境変数を設定する必要があるそう.
デスクトップ上にCygwinのホームディレクトリを置きたい場合
環境変数名: HOME
パス: /cygdrive/c/Users/username/Desktop
のようにWindows環境変数を設定する.