Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2018-05-29

MikuHatsune2018-05-29

ギブスサンプリング

Metropolis-Hastings サンプリングをやったので、ギブスサンプリングをやってみる。

ある変数X=¥{x_1,x_2,¥dots,x_{i-1},x_i,x_{i+1},¥dots,x_m¥} について、i 番目を取り除いたX_{-i}

X_i|X_{-i} ¥sim x_1,x_2,¥dots,x_{i-1},~~~x_{i+1},¥dots,x_mi 番目が抜けている)

で順次サンプリングして、その値を入れなおしてまたサンプリングする、を1¥dots m について行う。

 

結局、あるひとつi 以外のすべてを固定して、ひとつサンプリングすることになる。多次元正規分布の場合は、こちらを参考に

¥mu_{i|-i}¥leftarrow ¥mu_{i}+¥Sigma_{-i-i}^{-1}(X_{-i}-¥mu_{-i})

¥Sigma_{i|-i}¥leftarrow ¥Sigma_{ii}-¥Sigma_{-i-i}^{-1}¥Sigma_{-ii}

とする。ただし¥underbrace{¥Sigma_{-i-i}^{-1}}_{1¥times (m-1)}=¥underbrace{¥Sigma-ii}_{1¥times(m-1)}¥underbrace{¥Sigma_{-i-i}}_{(m-1)¥times(m-1)} である。なので¥mu_{i|-i}¥Sigma_{i|-i} はそれぞれ1¥times 1 の大きさというかスカラーで、各X_i はひとつずつサンプリングされる。

スクリプトここからパクる。

 

初期値が真の分布からは遠いところにしておいて、収束するまでにちょっと時間がかかるようにする。

Metropolis-Hastings(MH)では、真の分布に近づくまでに少し時間がかかる。サンプリングされた分布は、xy 同時にサンプリングしているので、移動方向は斜めである。

Gibbs (GB)は、各X_iサンプリングしているので、動き方は各x, y 方向にジグザグである。詳細釣り合い条件や移動確率が1であることはここでは取り上げないが、動きやすいのでMH より収束に向かいやすい。

この2次元の例では、一発目から真の分布に取り込まれている様子がわかる。

f:id:MikuHatsune:20180529171307p:image

 

多次元に拡張するが、図示することを考えて3次元にしてみる。MH もGB もひだり下のマゼンタのところからサンプリングを開始したが、HM は最初の50点の時点ではまだ収束していないが、GB は7点目くらいからもう収束している。

f:id:MikuHatsune:20180529171355p:image

MH の移動方向は斜めだが、GB では各点がサンプリングされるまでに、xyz でそれぞれX_iサンプリングしているので、カクカクになっている。

# 2D
# 初期値
x0 <- c(-5, 2)
iter <- 1500 # 繰り返し回数
x <- matrix(0, iter, 2) # ランダムウォークの座標
dratio <- rep(0, iter) # 確率密度比
out <- rep(0, iter) # 確率密度比にしたがって確率的に動いたか、動かなかったかの記録
x[1,] <- x0

# 2次元正規分布の相関
rho <- 0.7
sig <- matrix(c(1, rho, rho, 1), 2)
mu <- c(4, 5) # 2次元正規分布の各パラメータの平均

# Metrololis-Hastings
for(i in 2:iter){
  walk <- runif(2, -0.5, 0.5) # 歩く量
  v <- x[i-1,] + walk # 次の候補点
  dratio[i-1] <- dmvnorm(v, mean=mu, sigma=sig)/dmvnorm(x[i-1,], mean=mu, sigma=sig) # 確率密度比
  out[i-1] <- rbinom(1, 1, min(1, dratio[i-1])) # 動いたか、動かなかったか
  x[i,] <- x[i-1,] + out[i-1]*walk
}

kd1 <- kde2d(x[,1], x[,2], c(bandwidth.nrd(x[,1]), bandwidth.nrd(x[,2])), n=1000)
cols <- jet.colors(iter)
plot(x, type="p", pch=16, cex=0.6, col=2, xlab="x1", ylab="x2", main="Metropolis-Hastings sampling", xlim=xl, ylim=yl)
for(i in 2:iter){
  segments(x[i-1,1], x[i-1,2], x[i,1], x[i,2], col=cols[i])
}
points(x[1,1], x[1,2], pch="★", col=6)
points(mu[1], mu[2], pch="★", col=5)
contour(kd1, add=TRUE, col=1)

# Gibbs
X <- matrix(x0, nr=1)
for(j in 2:iter){
  x <- X[j-1,]
  for(i in seq(mu)){
    s <- sig[-i, i] %*% solve(sig[-i, -i])   # Σ_ab Σ_bb ^ -1
    # (PRML 2.81) μ_a|b = μ_a + Σ_ab Σ_bb ^ -1 (x_b - μ_b)
    mu_a_b <- mu[i] + s %*% (x[-i] - mu[-i])
    # (PRML 2.82) Σ_a|b = Σ_aa - Σ_ab Σ_bb ^ -1 Σ_ba
    sig_a_b <- sig[i, i] - s %*% sig[i, -i]
   # [Gibbs] x_a 〜 p(x_a|x_{-a}) = N(μ_a|b, Σ_a|b)
    x[i] <- rnorm(1, mu_a_b, sqrt(sig_a_b))
  }
  X <- rbind(X, x)
}
colnames(X) <- sprintf("V%d", seq(ncol(X)))

kd2 <- kde2d(X[,1], X[,2], c(bandwidth.nrd(X[,1]), bandwidth.nrd(X[,2])), n=1000)
cols <- jet.colors(iter)
plot(X, type="p", pch=16, cex=0.6, col=2, xlab="x1", ylab="x2", main="Gibbs sampling", xlim=xl, ylim=yl)
#lines(x, col=cols)
for(i in 2:iter){
  for(j in seq(ncol(X))){
    y <- rbind(replace(X[i-1,], 0:(j-1), X[i, 0:(j-1)]), replace(X[i-1,], 1:j, X[i, 1:j]))
    segments(y[1,1], y[1,2], y[2,1], y[2,2], col=cols[i])
  }
}
points(X[1,1], X[1,2], pch="★", col=6)
points(mu[1], mu[2], pch="★", col=5)
contour(kd2, add=TRUE, col=1)
# 3D
# 初期値
x0 <- c(-5, 2, -4)

iter <- 1500 # 繰り返し回数
x <- matrix(0, iter, length(x0)) # ランダムウォークの座標
dratio <- rep(0, iter) # 確率密度比
out <- rep(0, iter) # 確率密度比にしたがって確率的に動いたか、動かなかったかの記録
x[1,] <- x0
# 3次元の多次元正規分布を適当に作る
rho <- c(0.7, 0.6, 0.9)
sig <- diag(0, length(x0))
sig[lower.tri(sig)] <- rho
sig <- sig + t(sig)
diag(sig) <- rep(1, length(x0))
mu <- c(4, 5, 3)

# Metropolis-Hastings
for(i in 2:iter){
  walk <- runif(length(x0), -0.5, 0.5) # 歩く量
  v <- x[i-1,] + walk # 次の候補点
  dratio[i-1] <- dmvnorm(v, mean=mu, sigma=sig)/dmvnorm(x[i-1,], mean=mu, sigma=sig) # 確率密度比
  out[i-1] <- rbinom(1, 1, min(1, dratio[i-1])) # 動いたか、動かなかったか
  x[i,] <- x[i-1,] + out[i-1]*walk
}
colnames(x) <- sprintf("V%d", seq(ncol(x)))

cols <- jet.colors(iter)
plot3d(x, type="p", pch=16, size=0.01)
points3d(head(x, 50), size=5)
title3d("Metropolis-Hastings sampling", line=3)
for(i in 2:iter){
  segments3d(x[(i-1):i,], col=cols[i])
}
spheres3d(x[1,], col=6, radius=0.5)

# Gibbs
X <- matrix(x0, nr=1)
for(j in 2:iter){
  x <- X[j-1,]
  for(i in seq(mu)){
    s <- sig[-i, i] %*% solve(sig[-i, -i])   # Σ_ab Σ_bb ^ -1
    # (PRML 2.81) μ_a|b = μ_a + Σ_ab Σ_bb ^ -1 (x_b - μ_b)
    mu_a_b <- mu[i] + s %*% (x[-i] - mu[-i])
    # (PRML 2.82) Σ_a|b = Σ_aa - Σ_ab Σ_bb ^ -1 Σ_ba
    sig_a_b <- sig[i, i] - s %*% sig[i, -i]
   # [Gibbs] x_a 〜 p(x_a|x_{-a}) = N(μ_a|b, Σ_a|b)
    x[i] <- rnorm(1, mu_a_b, sqrt(sig_a_b))
  }
  X <- rbind(X, x)
}
colnames(X) <- sprintf("V%d", seq(ncol(X)))
# ギザギザの描画
Y <- X[1, , drop=FALSE]
for(i in 2:nrow(X)){
  for(j in 1:ncol(X)){
    y <- rbind(replace(X[i-1,], 0:(j-1), X[i, 0:(j-1)]), replace(X[i-1,], 1:j, X[i, 1:j]))
    Y <- rbind(Y, y)
  }
}

cols <- jet.colors(iter)
plot3d(X, size=0.01)
points3d(head(X, 50), size=5)
title3d("Gibbs sampling", line=3)
lines3d(Y, col=rep(cols, each=3))
spheres3d(x[1,], col=6, radius=0.5)

2018-05-25

ものすごいわかりやすかったKullback-Leibler divergence

情報量を-logQ(x) と置くと、

単調に減少する

ふたつの分布が独立ならば、加算的に扱える

というので、-logQ(x) と書いておく。

 

符号化してその平均長を考えると

¥sum-Q(x)logQ(x)エントロピーである。

さてここで、Q(x) はわからないけど、理論的にこれならば、上のエントロピーが成り立つ。しかし、何度も言うがQ(x) はわからないので、ひとまずP(x) を考える。これの平均長は

-¥sum xlogP(x)

だが、x¥sim Q(x) (実際にデータとして得ているx は、実態はわからないけれどもQ(x) からサンプリングされている)と書けるので、結局-¥sum Q(x)logP(x) である。

このふたつの、(真の分布Q(x) のときの平均長)-(真の分布がわからないけれどもとりあえずP(x) とおいた平均長)を計算すると

D(Q, P)=-¥sum Q(x)log¥frac{P(x)}{Q(x)}

となる。

これは常に0 以上で、{P(x)}{Q(x)}=X とすれば、y=X-1y=logX のふたつの関数を考えれば成り立ち、0 になるときは恒等的にP(x)=Q(x) のときである。

2018-04-25

MikuHatsune2018-04-25

ランダムウォークで乱数を生成する

機械学習のTA なのに機械学習素人なので機械学習統計の講義を聞いている。

多次元な確率分布から乱数をいい感じに取ってくる方法に、ランダムウォークを使う。

多次元な確率分布としては、単純に2次元正規分布、とする。二次元の正規分布は、x1 とx2 のそれぞれの変数が平均¥mu_1, ¥mu_2 で、変数間の関係として分散共分散行列¥Sigma=¥begin{pmatrix}¥sigma_1^2 & ¥sigma_1¥sigma_2 ¥¥ ¥sigma_2¥sigma_1 & ¥sigma_2^2¥end{pmatrix} である。

R ならrvtnorm パッケージが使える。

 

アルゴリズムとしては、二次元正規分布の確率密度関数g(x) として、ある適当な初期値x_0 から開始する。

あるx のとき、どれくらいランダムウォークするかは、一様分布U からランダムに決める。ここでは歩く量w¥sim U(-0.5, 0.5) とする。つまり、x方向もy方向も -0.5 から0.5 までの一様分布からサンプリングされる。

さて、ここでランダムサンプリングするかしないかは、以下のルールで決める。

あるx_{t} のとき一様分布からランダムウォークする量が決まって、v_{t} に移動したとすると、次の位置x_{t+1} は、x_{t+1}¥leftarrow v_{t} として位置を変えるか、x_{t+1}¥leftarrow x_{t} として動かないか、である。

この動く/動かない、を決めるのは、確率密度比¥frac{g(v_t)}{g(x_t)} より、a_t=¥textrm{min}¥{1, ¥frac{g(v)}{g(x)}¥}として確率を決めて、

P(x_{t+1}=v_t)=a_t:動く

P(x_{t+1}=x_t)=1-a_t:動かない

とする。

確率密度比¥frac{g(v_t)}{g(x_t)} が1 より大きいということは、g(v_t) のほうがg(x_t) よりもっともらしい、ということで、動くべき、ということになる。逆のときは動かないほうがいい、ということになるが、それでもちょっとは動くことも可能性としては残したほうがよい、という気分がある。というのが、(確率密度の)比は0 以下にはならないからである。

 

という感じで、適当な初期値、確率密度関数(ここではdrvnorm で得られる)、繰り返し回数、を指定すると、乱数ができる。繰り返しは10000回行った。

 

できた。

f:id:MikuHatsune:20180425202544p:image

 

さて、これは初期値依存なので、初期値がうまく設定されないと、いい感じにならない。ここで、初期値を少し離して実行してみる。

 

できた。

f:id:MikuHatsune:20180425202542p:image

 

ここで試している確率密度は単純なもので、初期値が目的の関数から離れていても、サンプリング回数が多くなればそこそこたどり着く、ことが多い。ここで、ある程度離れた初期値で、サンプリング回数も減らしてみる(1500回)。

 

最初の100サンプリングくらいでもう目的の関数部分にたどり着いているので、それなりにいい感じの乱数ができているが、上の図とはサンプル数が段違いなので、カーネル密度推定したものの楕円形(これは二次元正規分布で設定した¥Sigma の傾き具合に相当する)がグニャグニャである。

f:id:MikuHatsune:20180425202538p:image

 

上の図では、ランダムウォークの過程を各t のときで色付けしておいた。ここで、確率密度比¥frac{g(v_t)}{g(x_t)} (の対数)と、確率密度比の確率を用いて点が移動したのか否かの記録をプロットしている。

確率密度比¥frac{g(v_t)}{g(x_t)} (の対数)が非常に大きいもしくは小さいときは、最初の50点くらいで、それ以降はほとんど1、つまり¥frac{g(v_t)}{g(x_t)} がほぼ等しい状態になっている。この時点では、ランダムウォークの点は二次元正規分布に近いところをウロウロしている。ランダムウォークの移動量は、U(-0.5, 0.5) のけっこう狭いグリッド内なので、一発のランダムウォークで二次元正規分布から逸脱することはない。

¥frac{g(v_t)}{g(x_t)} が1 以上のときは、確率1 で移動することになるので、移動したことを示すマゼンタが記録されている。¥frac{g(v_t)}{g(x_t)} が1 未満のときは、確率¥frac{g(v_t)}{g(x_t)}で移動して確率1-¥frac{g(v_t)}{g(x_t)} でその場にとどまるので、¥frac{g(v_t)}{g(x_t)}=1(対数では0)より下の領域では、2色出ている。

f:id:MikuHatsune:20180425202536p:image

 

こういうのがRstan でのいわゆるwarm up や収束になる。

 

library(mvtnorm)
library(MASS)
x0 <- c(0.1, 0.2)
x0 <- c(-1, 2)

iter <- 1500 # 繰り返し回数
x <- matrix(0, iter, 2) # ランダムウォークの座標
dratio <- rep(0, iter) # 確率密度比
out <- rep(0, iter) # 確率密度比にしたがって確率的に動いたか、動かなかったかの記録
x[1,] <- x0

rho <- 0.7
sig <- matrix(c(1, rho, rho, 1), 2)
mu <- c(4, 5)

for(i in 2:iter){
  walk <- runif(2, -0.5, 0.5) # 歩く量
  v <- x[i-1,] + walk # 次の候補点
  dratio[i-1] <- dmvnorm(v, mean=mu, sigma=sig)/dmvnorm(x[i-1,], mean=mu, sigma=sig) # 確率密度比
  out[i-1] <- rbinom(1, 1, min(1, dratio[i-1])) # 動いたか、動かなかったか
  x[i,] <- x[i-1,] + out[i-1]*walk
}

kd <- kde2d(x[,1], x[,2], c(bandwidth.nrd(x[,1]), bandwidth.nrd(x[,2])), n=1000)

cols <- jet.colors(iter)
plot(x, type="p", pch=16, cex=0.6, col=2, xlab="x1", ylab="x2")
#lines(x, col=cols)
for(i in 2:iter){
  segments(x[i-1,1], x[i-1,2], x[i,1], x[i,2], col=cols[i])
}
points(x[1,1], x[1,2], pch="★", col=6)
points(mu[1], mu[2], pch="★", col=5)
contour(kd, add=TRUE, col=6)

pcols <- c("black", "magenta")
plot(dratio, log="y", type="n", xlab="sampling number", ylab="log likelihood")
for(i in 2:iter){
  segments(i-1, dratio[i-1], i, dratio[i], col=cols[i])
}
abline(h=1, lty=3)
points(seq(iter), dratio, cex=0.4, col=pcols[out+1], pch=16)
legend("topright", legend=c("動かなかったとき", "動いたとき"), pch=16, col=pcols, cex=2)

2017-11-20

MikuHatsune2017-11-20

logit の代わりにtanh を使ってはいけないのですか?

講義の最中にこんな質問があって、「logit は0-1 の範囲にあって、確率として扱いやすいから」という回答だったが、あとで

スケールが変わっただけで本質的には同じということがアナウンスされていた。

 

logit は

f(x)=¥frac{e^x}{1+e^x}

tanh は

g(x)=¥frac{e^x-e^{-x}}{e^x+e^{-x}}

となり、これらの関係は

g(x)=2f(2x)-1

となる。f(x)=g(x) となるのは、x_0=¥log(1+¥sqrt{2}) のときである。

 

logit が確率として頻用されるのは[0, 1] で収まるからだし、tanh が[-1, 1] で収まるのはdeep learning の活性化関数で頻用される。

また、こういう関数たちは連続なので微分が可能だから扱いやすいのだろう、と数学素人的に思う。

x <- seq(-5, 5, length=3000)
y <- cbind(logistic=gtools::inv.logit(x), tanh=tanh(x))
x0 <- log(1+sqrt(2))
cols <- c("blue", "red")
matplot(x, y, type="l", col=cols, lwd=5, lty=1)
legend("bottomright", legend=colnames(y), col=cols, lwd=5, cex=2)
points(x0, tanh(x0), pch=16, cex=2)

f:id:MikuHatsune:20171120191439p:image

2017-11-07

10万枚の胸部X線画像をCNNした

読んだ。

ChestX-ray8: Hospital-scale Chest X-ray Database and Benchmarks on Weakly-Supervised Classification and Localization of Common Thorax Diseases. IEEE CVPR 2017

3万人の患者から10万枚程度の胸部X線画像を入手し、CNNによりAtelectasis, Cardiomegaly, Effusion, Infiltration, Mass, Nodule, Pneumonia, Pneumothorax, Normal の8つの病気と正常の区別を行う。

F値ベースだと0.8-0.9 くらい、ただしNodule やMass だと0.5 とかになる。

AUCベースだと0.7-0.8 くらい、これもMass だと0.56 でやる意味あるの?ってなる。

判別器にはAlexNet, GoogLeNet, VGGNet16, ResNet50 を試していて、ResNet50 が総じてよさそうな性能。

 

Cardiomegaly 心肥大(AUC=0.8141) やPneumothorax 気胸(AUC=0.7891) はAUCが高く判別性能が高かったが、これらはぶっちゃけ人の目で見てもわかる。気胸はわからない場合もあるが、わからないような気胸は生命の危機が差し迫っているわけではないので翌日の受診でもよい。死にそうな気胸は息が止まるのでCNN に頼る必要がそもそも時間的にもない。

ではAI頼りのがんの検出のようなタスクはどうか。これはMass 大きめの腫瘤 (AUC=0.5609) Nodule もう少し小さめの腫瘤は検出性能が低かった。これは症例ごとのバリエーションが多いだろうからと言われている。

Pneumonia 肺炎(AUC=0.6333) も低い性能だったのは、データセット内に1% 程度しかなかったから、と言っている。

 

今回のデータセットは10万枚のうち8.4万枚が正常、ということで、なんらかの異常のある事前確率は15%程度になる。では、一般に健診でこんなに異常例があるかというと、たぶんない。疑いとして陽性をつけるかもしれないが、実際に病気である人は絶対に15%もいない。

この理由で、一般に気胸なんかは健診で見つけない(見つからない)。なんらかの症状があって、事前確率が高い状況で検査は行われる。

AIに期待されるがんの検出については、Mass やNodule の検出がたいした性能でなかったので、自分ならばこの程度の(診断責任を何も負わない)性能のAIに診断支援されたくないなぁというのが本音。

 

今回はF値とAUC とどちらも評価していたが、情報科学ではF値しか使われない印象がある。F値にはprecision とrecall が使われるが、これらは陽性に関わる項目(TP かFP かFN)しか使われない。

本当の理由はよくわからないが、ネットでググってよく出てくるのが、例えば書類を検索するときに、検索した書類が目的のものか(precision)、存在するすべての書類のうちどれくらい検索できたか(recall)という、書類を見つけ出すことにしか興味がなさそうである。

一方で医学では、病気を病気と診断すること(感度)と、病気でないことを病気でないと診断すること(特異度)が必要になってくる。

F値だけ出されると0.8-0.9程度といい性能だと勘違いしそうだが、8割以上が陰性例なので(というか臨床検査なので)AUC でも評価しているのは好感がもてた。性能についてはまったく期待できないが。

 

結局のところ、Net組み上げました選手権になるか、超高性能GPU使って入力画像増やしまくりました選手権になるだけだと思う。はやくAIが仕事奪いに来いよ()

リンクから生データが入手できるので、CNN もできるし画像読影の勉強もできる。