ひらめの日常

日常のメモをつらつらと

RでKLダイバージェンスを描画する

今回はRでKLダイバージェンスを描画します。
KLダイバージェンスとは、カルバック・ライブラー情報量の略称。2つの確率分布の差異を測る距離的尺度として捉えることができます。

と言われても実感がわかないと思うので、実際に描画して感覚を掴んでいきたいと思います。
実装はこちらにあります。
7/5:リンク先更新しました。

github.com

定義

まずはエントロピーの説明をし、その上にあるKLダイバージェンスの定義を述べます。

エントロピー

連続関数のエントロピーは、情報量の期待値という定義なので次のように表せます。

\displaystyle
H(X) = -\int_{-\infty}^{\infty}p(x)\ln{p(x)} dx
直感的には、次に得られる情報量がどれだけの価値がある可能性があるかを示している数値になります。
(情報量とは、確率的に低いものを得た時に情報量を多く得たと判断するものです。)

KLダイバージェンス

今、正しい確率分布  p(x) と、その確率分布を近似したモデルである  q(x) があるとします。
この時に、 q(x) によって真の値を推測するために追加で必要な情報量の平均として定義されます。
具体的には、次のように表されます。「 q(x)エントロピー-  p(x)エントロピー」と捉えるとわかりやすいかもしれません。
{\displaystyle
KL(p||q) = -\int_{-\infty}^{\infty}p(x)\ln{\frac{q(x)}{p(x)}}dx
}

描画してみる

対象とするモデル

以下のように、平均mean2=2, 標準偏差sd2=2の正規分布を対象とします。

X = seq(-10, 10, length = n)
Y4 = dnorm(X, mean = mean2, sd = sd2)
plot(X,
Y4,
type = "l",
col = "blue")

f:id:thescript1210:20180501005147p:plain

KLダイバージェンスの図示

左の図は平均値を変化させた時のKLダイバージェンスの値。平均値が対象モデルと同じ2の時にKLダイバージェンスは最小値を取ることがわかります。
右の図は標準偏差を変化させた時のKLダイバージェンスの値。標準偏差が対象モデルと同じ2の時にKLダイバージェンスは最小値を取ることがわかります。

これらのことから、KLダイバージェンスは、モデル間のパラメータを用いて、モデル間の距離を表していると捉えることができますね。

f_kl_divergence = function(m1, m2, sd1, sd2) {
  log(sd1 / sd2) + (sd2 ^ 2 + (m1 - m2) ^ 2) / (2 * sd1 ^ 2) - 1 / 2
}
y_mean = c()
x = seq(-4, 4, length = 100)
for (i in x) {
y_mean = append(y_mean, f_kl_divergence(i, mean2, sd2, sd2))
}

y_sd = c()
x2 = seq(1, 3, length = 100)
for (i in x2) {
   y_sd = append(y_sd, f_kl_divergence(mean2, mean2, i, sd2))
}

layout(matrix(1:2, ncol=2))
plot(x, y_mean, type="l")
plot(x2, y_sd, type="l")

f:id:thescript1210:20180501005336p:plain

(番外編)Rのライブラリを用いてみる

Rには、KLダイバージェンスにしたがって乱数を生成し、元モデルのパラメータを推定する関数があります。
その関数であるKL.divergence()を用いて乱数を発生させ、それをプロットしました。

するとその結果は先ほどplotした理論値に近いグラフを描画することができました。

library(FNN)
set.seed(100)

kl_mean = c()
Y1 =  rnorm(n * 10, mean = mean2, sd = sd2)
for (i in seq(-4, 4, length = 100)) {
  Y2 = rnorm(n * 10, mean = i, sd = sd2)
  divs = KL.divergence(Y1, Y2, k = 2)
  kl_mean = append(kl_mean, mean(divs))
}

kl_sd = c()
for (i in seq(1, 3, length = 100)) {
  Y3 = rnorm(n * 10, mean = mean2, sd = i)
  divs2 = KL.divergence(Y1, Y3, k = 2)
  kl_sd = append(kl_sd, mean(divs2))
}

layout(matrix(1:2, ncol=2))
plot(seq(-4, 4, length = 100), kl_mean, type="l")
plot(seq(1, 3, length = 100), kl_sd, type="l")

f:id:thescript1210:20180501005418p:plain

Rで2次元正規分布を描画する

自分の実装はgithubにあります。

github.com

D次元正規分布

1次元正規分布はよく統計学とかでも出て来ますが、今回は2次元正規分布確率密度関数を描画したいと思います。

定義

{ \displaystyle
f(x) = \frac{1}{(2\pi)^{\frac{D}{2}}|\Sigma|^{\frac{1}{2}}}\exp(-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu))
}

{D}:次元
{x}:D次元ベクトル
{\Sigma}:xの分散共分散行列
{\mu}:D次元の分散ベクトル

以後D=2の場合を扱っていきます。

散布図

とりあえず確率密度関数を描画する前に、2次元正規分布に従う散布図をプロットしてみます。

library(mvtnorm)
sigma = matrix(c(2, 0.8, 0.8, 2), ncol = 2)
rand = rmvnorm(n = 1000, mean = c(0, 0), sigma)
x1 = rand[, 1]
x2 = rand[, 2]
plot(x1, x2)

mvtnormが多次元正規分布を描画するときには必要になります。
流れとしては、

  1. sigmaという分散共分散行列の定義
  2. rmvnormの引数に、平均と分散共分散行列を渡す
  3. rmvnormは渡された次元数を元に、多次元正規分布に従った乱数を生成

f:id:thescript1210:20180425000159p:plain

確率密度関数

それでは、確率密度関数{f(x)}を描画していきます。

x1 = seq(-3, 3, length = 50)
x2 = x1
f = function(x1, x2) {
  dmvnorm(matrix(c(x1, x2), ncol = 2), mean = c(0, 0), sigma = sigma)
}
z = outer(x1, x2, f)
z[is.na(z)] = 1
op = par(bg = "white")
persp(x1, x2, z, theta = 30, phi = 30, expand = 0.5, col = "orange")  

今回は、dmvnormという関数を用います。

  1. x1, x2に軸用の値を代入
  2. dmvnormの第一引数に軸ラベル、そして平均と分散共分散行列を渡す
  3. dmvnormは渡された次元数を元に、多次元正規分布を生成
  4. 多次元正規分布の出力をzに代入しperspでグラフ描画

こんな感じになります。
f:id:thescript1210:20180425000954p:plain

参考

こちらのサイトなどを参考にさせていただきました。
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html

Shift-JISのエクセルファイルの文字化けを直す

何回も同じことを調べていて覚えないので備忘録的に残します。
ファイルの文字コードを調べる

>> nkf -g hoge.csv
Shift-JIS

ファイルの文字コードをUTF-BOMに変更する
UTFの文字コードにすると、エクセルは適切に読み込んでくれないので、UTF-BOMに変換する必要があるようです。

>> nkf --overwrite --oc=UTF-8-BOM hoge.csv
>> nkf -g hoge.csv
UTF-8

メモリの種類を10分でまとめる(ROM, RAM)

◎ROM

Read Only Mmory
その名の通り基本動作は読み出しのみ。
次のようにROMには様々な種類があります。

○マスクROM

メモリLSI(データを記憶に特化した回路のこと)の設計時に内容が決まる。ユーザーによる消去・書き込みは不可能。

○PROM ROM

Programmable ROM
ユーザーによる消去・書き込みが可能。

ヒューズROM

ユーザーが一度だけ書き込みが可能。以後の消去・書き込みは不可能。

EPROM

Erasable PROM
ユーザーによる消去・書き込みが何度でも可能。

  • UVEPROM(紫外線で消去・書き込み)
  • EEPROM(電気で消去・書き込み)
  • フラッシュメモリ(電気で消去・書き込み。ブロック単位で高速に処理可能)

◎RAM

Random Access Mmory
基本動作は読み出しと書き込み。
次のようにRAMにも様々な種類があります。

SRAM

Static RAM
セルはフリップフロップで、リフレッシュ(後述)不要。
高速で、レジスタファイルやキャッシュメモリなどに用いられる。

DRAM

Dynamic RAM
セルはコンデンサ
コンデンサなので、放置しておくと電荷が流れてしまう。そこで、ある時間ごとに電荷の状態を確認して再び電荷を加え直したりしなければならない。このことをリフレッシュという。

高速モードDRAM

連続する列アクセスを高速化。

SDRAM

Static DRAM
クロック同期とアクセスのオーバーラップによる高速化。

RDRAM

Rambus DRAM
データ幅の縮小と直列接続による高速化。


こちらのpdfファイルを参考にしました
http://www.mtl.t.u-tokyo.ac.jp/~sakai/hard/hard2.pdf

コマンドラインからpickleファイルの中身確認

Pythonファイルでは

pythonではpickle.dump()で書き込み、pickle.load()で読み込みます。

# 書き込み
with open("hoge.pkl", mode="wb") as f:
    pickle.dump({"hoge":"fuga"}, f)

# 読み込み
with open("hoge.pkl", mode="rb") as f:
    hoge = pickle.load(f)

コマンドラインでは

しかし、これでは手軽にpickleファイルの中身を確認することができません。
pythonファイルを別に作ってそれを毎回指定して実行するのも面倒な気がします。

すると、コマンドラインからでも次のようにすれば読み込めるとのこと。

$ python -m pickle hoge.pkl
{hoge: fuga}

これだけ。

こちらの記事を参考にさせていただきました。
qiita.com

python+scikit-learnで多項式曲線近似をリッジ回帰で求める

はじめに

実装はこちらに全て載っています。
7/5:リンク先変更しました。 github.com

多項式曲線近似とは

多項式曲線近似とは、関数があったとき、それを次のような多項式で近似することです。
{ \displaystyle
y(x, w) = w_0 + w_1x + w_2x^{2} + ... + w_Mx^{M} = \sum_{j=0}^{M} w_jx^{j}
}
ここで、Mはパラメータで、自分で近似する多項式の次数を決定します。一般的に次数が高すぎると過学習してしまい、学習データセットに対する精度は上がりますが、汎化性能については低くなってしまうことが多いです。

これを一般的に、線形モデルと言います。

係数の求め方

ここで、線形モデルの各係数wを求める必要があります。
多項式近似して予測した値と、実際の値との誤差関数を最小化するという方針で求めます。

一般的な回帰モデル

誤差関数は以下が用いられることが多いです。二乗誤差の最小値を求めに行くので、最小二乗法と呼ばれます。

{ \displaystyle
E(w) = \frac{1}{2} \sum_{n=1}^{N} \{y(x_n, w) - t_n\}^{2}
}
ここで、y(x_n, w)多項式によって予測された推定値。t_nは正解の値です。この誤差関数をw_j微分して最小値を求めます。

正則化した回帰モデル

一般的な回帰モデルだと、元の関数から大きくずれることがあります。それは、学習データセットに対して適合し過ぎてしまい、大きい係数が出てきてしまうためです。

そこでリッジ回帰モデルでは、罰則項を付加することにより過学習を防ぎます。この作業は正則化(regularization)と呼ばれています。

具体的には以下のような誤差関数によって回帰モデルを求めます。
最後に係数を二乗で足しているのがわかりますね。
{ \displaystyle
E(w) = \frac{1}{2} \sum_{n=1}^{N} \{y(x_n, w) - t_n\}^{2} + a|w|^{2}
}
ここで、aはパラメータで自分で値を設定します。

実装

元の関数

python+sklearnで実装しました。
元の関数、そしてそこから外れた値を生成するためにrandomな値を加算する関数を定義します。

def function(x):
    return 5 * np.sin(x * np.pi / 10)

def function_with_random(x):
    return function(x) + np.random.uniform(-5.0, 5.0)

多項式の定義

sklearnには、PolynomialFeaturesというクラスがあり、今回6次の多項式を定義しました。

from sklearn.preprocessing import PolynomialFeatures
polynomial_features = PolynomialFeatures(degree=6)
X_poly = polynomial_features.fit_transform(X)

このfit_transform()ですが、これは定義した多項式に対して、引数で与えられたXの値を代入した配列が入っています。今回は以下のようになります。

[[1, 0, 0, 0, 0, 0, 0],  # x = 0
[1, 1, 1, 1, 1, 1, 1, 1],  # x = 1
[1, 2, 4, 8, 16, 32, 64],  # x = 2
...
]

学習

次にやることは、先ほど各 1, x, x^{2}, x^{3}...の値を計算したので、それらの線形結合で予測値を出すということです。

from sklearn.linear_model import Ridge
poly_reg = Ridge(alpha=10, fit_intercept=False)
reg  = poly_reg.fit(X_poly, y_train)  # X_polyの線形結合によりyの値を学習する

予測

seabornを使ってグラフを描画しました。

y_ridge = poly_reg.predict(X_poly)  # 学習したモデルを用いて予測する。
sns.pointplot(x=X.reshape(len(X)), y= y_ridge, color='R', markers="")

赤がリッジ回帰によって得られた曲線。青が比較のために実装した正則化しない線形回帰によって得られた曲線です。
確かに外れ値の値を受けにくいような気がします。
f:id:thescript1210:20180409010132p:plain
モデルの平均二乗誤差は次のようになり、確かにリッジ回帰モデルでは過学習が抑制され、汎化性能を保っていることがわかります。

rigde_train_error: 3.8009207105204244
rigde_test_error: 2.665180445057409
linear_train_error: 1.987241639408757
linear_test_error: 3.6306109275834806

ChainerCVでSSDを学習させる時の注意点

はじめに

この記事は、自分がChainerCVを用いてSSDを学習させる時にハマった点や注意した点を紹介しています。
ChainerCVとは、その名の通りChainerをベースにコンピュータビジョンのタスクに対応させたものとなります。

SSDとは

簡単にSSDの説明をすると、SSDとはSingle Shot MultiBox Detectorの略で、

①VGGをbase networkとして異なる大きさのフィルターをかけていき
②異なるアスペクト比、解像度のbounding boxに対応した
③精度74.3%, 速度59FPSを記録した高い精度で高速

な物体検出方法です。それまでの先行研究に比べて学習時に物体らしさのを物体の分類と同時に学習するために、今までよりも高速な学習が可能になっています。

論文はこちら。
[1512.02325] SSD: Single Shot MultiBox Detector

また、こちらのスライドまとめがSSD含めた物体検出(object detection)の理解に非常に役立ちました。

www.slideshare.net

ChainerCVとは

ChainerCVは、ディープラーニング用ライブラリChainerのコンピュータービジョンに対応したもので、物体検出の手法が数多く実装されています。
GitHub - chainer/chainercv: ChainerCV: a Library for Deep Learning in Computer Vision

参考にした学習コードは、ChainerCVのexamplesにあるtrain.pyです。
github.com

注意点

MultiProcessIteratorのデッドロック

よし、gpuで学習を始めよう!と思うわけですが、学習が一向に進みません。SerialProcellIteratorを用いた時には正常に動くので、MultiProcessIterator特有の挙動のようです。

ChainerCVのissueを漁っていると同じような状況ハマっている方を見つけました。
github.com
どうやらMultiProcessIteratorが内部でcv2.resize()を呼んでおり、これが別スレッドを立てるためにデッドロックに陥るとのことです。

コメント通り、util.pycv2.resize()の直前にcv2.setNumThreads(0)を追加したところ、学習が開始されました。

Validationの適切なログ出力

ここでいう「適切」とは、自身に秘帖なタイミングで適切に出力するように変更しようということです。

デフォルトでは、10,000iterationごとにtestデータセットを用いてvalidation accuracyをlog出力しています。学習を120,000iteration回すので、それでも良いかもしれませんが、自分は頻繁にログを確認したかったので次のように変更しました。

trainer.extend(
        DetectionVOCEvaluator(
            test_iter, model, use_07_metric=True,
            label_names=voc_bbox_label_names_with_hands),
        trigger=(1000, 'iteration'))  # 1000iterationごとにvalidationのログ出力

頻繁にvalidation accuracyを出力する時ですが、かなり時間がかかるのでそこは留意したほうが良いかと思います。

使用できる学習済みモデルはimagenet

ChainerCVで実装されているSSDはVOCデータセットの20クラス分類に対応しています。
ですので、自作のデータセット含めた21クラス分類の転移学習(fine-tuning)がVOCデータセットの重みを用いてできるかと思ったのですが、現状はできません

こちらのissueにもあるように、20クラス分類以外の場合はVOCデータセットの重みは用いることができず、現在実装中のようです。
github.com

自作データセットを用いて20クラス分類以外の学習をさせたい時には、pretrained_modelはimagenetのままで学習させましょう。

if args.model == 'ssd300':
        model = SSD300(
            n_fg_class=21,  # 用途に応じて変更
            pretrained_model='imagenet')
    elif args.model == 'ssd512':
        model = SSD512(
            n_fg_class=21,  # 用途に応じて変更
            pretrained_model='imagenet')