『わけがわかる機械学習』中谷秀洋（著）の書評

僕が中谷さんと初めて会ったのはみどりぼんの読書会で、初めて話したのは岩波DSの打ち合わせだったと思います。今でもそんなに親しくはないと思います。しかし、中谷さんのブログは10年ぐらい前から読んでいました。自然言語処理を中心とする機械学習に関連する理論（の解釈）・論文レビュー・数値実験の記事が多く、他のブログでは見られない独特かつ理解の深い内容で、毎日勉強させてもらっていました。今でも何度も読むべきブログです。その中谷さんが機械学習についてまるごと一冊書いたものが本書になります。もともと買うつもりでしたが、献本いただいたので簡単にご紹介いたします。

わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する

目次は以下になります。

• 0章： はじめに
• 1章： 機械学習ことはじめ
• 2章： 確率
• 3章： 連続確率と正規分布
• 4章： 線形回帰
• 5章： ベイズ確率
• 6章： ベイズ線形回帰
• 7章： 分類問題
• 8章： 最適化
• 9章： モデル選択
• 10章： おわりに
• 付録A： 本書で用いる数学

特長

一読して最初の感想は前半（2章～3章, 5章）が「確率および分布の初中級者向けの完ぺきな入門」、後半（4章～7章）が「わけがわかるPRML」あるいは「声に出して読みたいPRML」です。

残りの部分は、深層学習を含む機械学習全般（ところどころ深層学習にも言及しています）で必須の内容（モデルとは、交差検証、など）を扱っています。

理論の本にありがちな「天下りで定義や式が与えられて、わけわからんけど式変形していく」ではなく、「経験から確率は〇〇という条件を満たしてほしいよね。仮に〇〇として議論を進めるとどうだろう？」という感じに、ボトムアップで話が展開されていくので納得しながら読みやすいです。

5章は4章を拡張するために導入され、8章は7章の問題を解くために導入されます。章のタイトルだけ見ると一見関係性が不明ですが、流れるように議論をすすめるためにこの配置になっています。

数式に関しては0.1節に

本書は「数式がなくてもわかる本」ではありません。

と書いてあるように大学学部レベルの数式が出てきます。しかし、遠慮なく式が出てくるのは4章以降ですし、式変形の内容や考え方が非常に丁寧に書いてあるので、じっくり読めば理解できて力になることは間違いありません。

以降では部分的に補足します。

2～3章

0, 1, 9, 10章は合計して25ページぐらいなのに、2章と3章は合計して65ページもあって力の入れようが分かります。一見普通の章タイトルですが、内容はめちゃくちゃいいです。特に例題をまじえて「同時分布・周辺分布・条件付き確率分布」を説明しているあたりがいいです。

僕のブログの以前の記事で

色々読んでみましたが、現在決定版と言えるものは存在しないように思えました。個人的には、シグマと積分の復習、場合の数・数え上げの方法、確率、確率変数、確率密度、度数分布とヒストグラム、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換、検定、散布図と箱ひげ図、回帰、相関あたりをRやPythonなどを使いながらシンプルに説明していく本があるといいと思うのですが、なかなかバランスのとれたいい本がありません。初歩の初歩しか説明してない、グラフが少ない、検定にページを割きすぎ、分厚い、ちょっと難しいなどの不満点があります。立ち読みして自分にあった本を選ぶのがいいと思います。

と書きましたが、「確率、確率変数、確率密度、代表値・平均・分散、確率分布、同時分布、周辺分布、確率変数の変数変換」に関して僕が想定する内容と比べて200%ぐらいの完成度です。これからは「そのあたりがわかんなかったら、この本の2,3章読んでおいてね」と言うことにします。

4～7章

4.1節 最小二乗法の最後の段落（p.88）で

ここまで、自然な流れに従って「正解」を出したように感じるかもしれません。しかしこの短い話の中で大きな仮定を4つも使っています。仮定1：関数の形が1次式。 仮定2：二乗和誤差が最小＝良い答え。仮定3：誤差独立の仮定。仮定4…（略）

というように前提条件をきっちり説明したうえで、「どこが拡張しやすいか」「拡張したらどうなるか」を続けて論じていきます。一貫してこのような理屈の説明が出てくるのが本書の特徴と思います。

7章は中谷さんの専門の自然言語処理の例題で話がすすむためか、モデルの考え方や苦労などもちらほらあって、勉強になるのはもちろんとても楽しく読める章になっています。

4～7章は、PRMLの3～4章を読もうとしたけど挫折した人、読んだけどイマイチ分かってない人に特にすすめたいです。

Enjoy!

NIMBLEでノンパラベイズを試したメモ

NIMBLEというRのライブラリがあります。BUGS言語風の文法でC++にコンパイルされるタイプの確率的プログラミング言語です。実装されている推定のアルゴリズムはここに書いてあります。MCMCの他にも以下のようなアルゴリズムがデフォルトで実装されており、実行速度もかなり速いです。

• particle filter（Sequential Monte Carlo, SMC）
• 空間モデルのCARモデル（Gaussian Markov Random Fieldsと等価）
• ノンパラベイズのChinese Restaurant Process（CRP）とStick-breaking Process（SBP）

ここでは、ノンパラベイズのCRPとSBPを試してみます。参考にしたのは公式ドキュメントです。

データは以下の以前の記事と同じで、モデルも似たようなものです。

statmodeling.hatenablog.com

はじめにCRP版、次にSBP版を紹介します。

Chinese Restaurant Process（CRP）

library(nimble)

code <- nimbleCode({
  group[1:N] ~ dCRP(alpha, size=N)
  alpha ~ dgamma(1, 1)
  for(c in 1:C) {
    mu_mix[c] ~ dnorm(20, 100)
    s2_mix[c] ~ dinvgamma(5, b)
  }
  b ~ dgamma(0.01, 0.01)
  for(n in 1:N) {
    Y[n] ~ dnorm(mu_mix[group[n]], var=s2_mix[group[n]])
  }
})

set.seed(1)
load('data-nonpara-bayes.RData')
constants <- list(N = data$n, C = data$C)
data <- list(Y = data$velocity)
inits <- list(mu_mix = seq(0, 20, len = constants$C),
              s2_mix = rep(4, constants$C),
              group = sample(1:10, size = constants$N, replace = TRUE),
              b = 1, alpha = 1)

fit <- nimbleMCMC(code = code, constants = constants,
                  data = data, inits = inits,
                  nchains = 4, niter = 10000, nburnin = 2000,
                  summary = TRUE, samplesAsCodaMCMC = TRUE,
                  monitors = c('group', 'mu_mix', 's2_mix' , 'alpha'))

library(coda)
pdf(file='fit-traceplot-CRP.pdf')
traceplot(fit$samples)
dev.off()

• 4行目の group[1:N] ~ dCRP(alpha, size=N) でCRPに従うサンプリングができます。専用のアルゴリズムになっているとのことで、BUGSやStanに比べると動作がだいぶ高速です。
• 6行目の C はクラスター数で、Nよりはだいぶ小さく、ただし、クラスターとクラスターがつぶれて合体しないように少し大きめにとる必要があります（nimbleを実行する際にその旨のwarningが表示されます）。データ（data-nonpara-bayes.RData）に含まれていますが、ここでは20に設定しています。
• 20行目　Stanと異なり、全てのパラメータの初期値をもれなくきちんと与えておかないと、尤度が無限やNAになってRごと落ちたりしますので注意。
• 28行目 samplesAsCodaMCMC = TRUEを設定しておくと{coda}ライブラリで扱いやすい形のサンプリング結果となり、31～34行目のようにしてtraceplotを出力することができます。

その他にも、コンパイルとサンプリングの実行を別にできたりしますが、ドキュメントが若干分かりにくいです。

Stick-breaking Process（SBP）

library(nimble)

code <- nimbleCode({
  for(c in 1:(C-1)) {
    v[c] ~ dbeta(1, alpha)
  }
  alpha ~ dgamma(1, 1)
  w[1:C] <- stick_breaking(v[1:(C-1)])
  for(c in 1:C) {
    mu_mix[c] ~ dnorm(20, 100)
    s2_mix[c] ~ dinvgamma(5, b)
  }
  b ~ dgamma(0.01, 0.01)
  for(n in 1:N) {
    group[n] ~ dcat(w[1:C])
    Y[n] ~ dnorm(mu_mix[group[n]], var=s2_mix[group[n]])
  }
})

(...CRP版と同じ...)

inits <- list(mu_mix = seq(0, 20, len = constants$C),
              s2_mix = rep(4, constants$C),
              group = sample(1:10, size = constants$N, replace = TRUE),
              v = rbeta(constants$C, 1, 1),
              b = 1, alpha = 1)

(...CRP版と同じ...)

CRP版とほとんど同じです。

どちらの場合においても、サンプリング結果からクラスター数や密度を算出する関数を自作する必要があります。またモデルのコードにバグが含まれていると、コンパイルエラーではなくRごと落ちたりするのは若いライブラリのご愛敬です。多少不便でもMCMCの性能はかなり良いと思うので積極的に使いたいと思います。

Enjoy!

Stanで1変数の積分を実行する

Stan2.19でモデルの内部で1変数の積分を実行する機能が追加されましたので簡単に紹介します。Stan2.19がCRANにないとお嘆きのアナタ！GitHubからインストールすればよろし。Windowsならば、Benさんがここに書いてくれた方法でバイナリを直接インストールすることもできます。ただし、他のパッケージとの依存関係は壊れてもしーらん。

積分する関数名は現在のところintegrate_1dです。使い方はStanのマニュアルに載っているかなと期待したけどまだ載ってませんでした。reference manualには少しだけ触れられているけど、これだけではさっぱり分からないです。こういう場合は、StanのdiscouseかGitHubのリポジトリで検索して、手あたりしだい見ていくと知見が蓄積されます。結局、GitHubのリポジトリ内で新しい関数のテストをしているディレクトリがあり、そこに使用例がありました。

この記事では階層モデルの場合に、積分を使って対数尤度を算出する例を扱います。階層モデルの場合には複数の予測したい量が考えられて、それに応じた対数尤度（および情報量規準）が考えられるのでした。

statmodeling.hatenablog.com

ここでは上記の記事でいうところの「モデル2の(3)」に対応する予測を扱います。すなわち、別の新しいクラスができて新しく1人が加わる場合です。この場合、予測分布はグループの平均を積分消去した分布になります。

先にデータ生成およびStanを実行するRコードから説明します。

library(rstan)

G      <- 10
N_by_G <- 2
N      <- N_by_G * G
GID    <- rep(1:G, N_by_G)

set.seed(2)
Mu <- rnorm(G, mean=0, sd=10.0)
Y  <- rnorm(N, mean=Mu[GID], sd=2.0)

data <- list(N=N, G=G, N_by_G=N_by_G, GID=GID, Y=Y)
stanmodel <- stan_model(file='model/model.stan')
fit <- sampling(stanmodel, data=data, iter=3500, warmup=1000, seed=123)

library(loo)
ms <- rstan::extract(fit)
waic  <- waic(ms$log_lik)$waic/(2*N)
looic <- loo(ms$log_lik)$looic/(2*N)

• 3行目： グループの数です。
• 4行目： 1グループに何人いるかです。
• 5行目： 全部で何人いるかです。ここでは20人です。
• 9行目： 各グループの平均を生成しています。
• 10行目： 9行目で生成した値を平均として、各グループに含まれる人たちの値を生成しています。
• 16～19行目： 対数尤度のサンプリング結果から(3)の場合のWAICとLOOICを算出しています。

Stanコードは以下です。

functions {
  real f(real x, real xc, real[] par, real[] y_r, int[] y_i) {
    return exp(normal_lpdf(x   | par[1], par[2]) +
               normal_lpdf(y_r | x,      par[3]));
  }

  real log_lik_marginal(real[] Y, real[] par) {
    return log(integrate_1d(f, par[1]-6*par[2], par[1]+6*par[2], par,
                            Y, {0}, 1e-3));
  }
}

data {
  int G;
  int N;
  vector[N] Y;
  int GID[N];
}

parameters {
  real mu0;
  real<lower=0> s0;
  vector[G] mu;
  real<lower=0> sigma;
}

model {
  mu0 ~ student_t(6, 0, 10);
  s0 ~ student_t(6, 0, 10);
  mu ~ normal(mu0, s0);
  Y ~ normal(mu[GID], sigma);
}

generated quantities {
  vector[N] log_lik;
  real mu_pred = normal_rng(mu0, s0);
  real y_pred = normal_rng(mu_pred, sigma);
  for (n in 1:N)
    log_lik[n] = log_lik_marginal({Y[n]}, {mu0, s0, sigma});
}

• 2～5行目： 被積分関数です。現状ではこの引数タイプ（real, real, real[], real[], int[]）に限るようです。1番目の引数は積分で動く変数、2番目は後述、3番目の引数はパラメータの配列、4番目の引数は実数値をとるデータの配列、5番目の引数は整数値をとるデータの配列です。2番目の引数はどんなときに使うべきか不明だけど、たぶん3番目の引数とは切り離して考えたいパラメータを入れるのに使うのかなぁ…。
• 7～10行目： 積分して対数尤度を求める関数の定義です。integrate_1d関数の、1番目の引数は関数名、2番目の引数は積分範囲の下限、3番目の引数は積分範囲の上限、4番目の引数は関数に渡すパラメータの配列、5番目の引数は関数に渡す実数値をとるデータの配列、6番目の引数は関数に渡す整数値をとるデータの配列、7番目の引数は積分の相対的な誤差をどれくらいに抑えるか、です。最近のStanでは{}で囲むと配列となり、[]で囲むとvectorとなります。
• 39行目： ここで対数尤度を求めています。

実際に実行してみると、Rでデータを作成するときの乱数をset.seed(1)の場合には情報量規準の算出までいきますが、set.seed(2)の場合には以下のエラーが出てサンプリングが止まってしまいます。

...(中略)...
Chain 2: 
[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                                                                      
[2] "  Exception: Exception: integrate: error estimate of integral above zero 1.80655e-006 exceeds the given relative tolerance times norm of integral above zero  (in 'model2f8834cb4d93_model' at line 8)"
[3] "  (in 'model2f8834cb4d93_model' at line 39)"                                                                                                                                                           
[1] "error occurred during calling the sampler; sampling not done"

Stanのintegrate_1dは内部的にはBoostのintegrate_1dを使っているようですが、どうもやや不安定な気がします。GSLのCQUADの方が安定しています（{RcppGSL}で使う例はこちら）。また、尤度を算出するような場合、どこかに対数をいれて積分しないと値が小さくなりすぎて数値的に不安定になります。以前のlog_sum_expとシンプソンの公式を使う方法は一例ですが、特化した効率の良い方法がありそうです。今後に期待します。いや、お前がやるんだよ。