やってみた。プロ棋士の強さ推定

階層ベイズモデルで勝敗データからプロ棋士の強さを推定する,StatModeling Memorandumという実験があります。大変面白いので私もやってみました。

データ収集

まず将棋のプロ棋士の対局結果を将棋連盟 棋士別成績一覧から収集しました。こちらのサイトは1964年度以降の棋士公式戦対局結果がまとまっています。米長邦雄の四段昇段が1963年、中原誠の四段昇段が1965年です。大昔ですね。素晴らしいですね。以下の方針でスクレイピングしました。

  1. データベース中のプロ棋士同士の対局をすべて収集する
  2. 対局者のどちらかあるいは両方が四段昇段前(奨励会員もしくはアマチュア)であっても収集対象とする
  3. 不戦勝、不戦敗は除外

2017年7月30日時点で、対局総数は103849局、棋士総数292人でした。

また、参考用にプロ棋士の生年月日および四段昇段日も収集しました。こちらは将棋連盟のデータベース将棋順位戦データベースWikipediaの個別ページを基にしました。将棋順位戦データベースも暗記したいほどの素晴らしいサイトです。

モデル

Memorandumの記事で使われたモデルのキモは次の2点だと思います。

  • 棋士のskill値+勝負ムラの大小で対局の勝敗をモデル化する
  • 棋士のskillの推移を自己相関モデルでモデル化する

モデルはBUGS文法で書かれており、とても複雑に見えます。BUGSはStan以上に使い方が分からないので、私はMemorandumのモデルをせめてStanに移植できないかと試行錯誤しましたが、何をどうしてもサンプリング結果が収束しませんでした。残念です。素人には高いハードルでした。

しょうがないのでモデルを変更しました。対局の勝敗は、勝者の勝率pをskill値の差(勝者-敗者)によるロジットモデルで表します。勝率pはそのまま尤度となります。勝負ムラは考慮しません。


logit(p) = a \times \Delta skill

このロジットモデルというのは今回初めて知りましたが1、いかにも便利で使い出がありそうですね。

skillの自己相関は下のように定めました。iは年度のインデックスです。sは全棋士全年度で共通としてます。また、これを正しく自己相関と呼んでいいのか、私は知りません。


skill_i \sim \mathcal{N}(skill_{i-1}, s^2)

aが1/400 * log(10)の場合、スキル差と勝率の関係がイロレーティングと同じになります。つまりこのモデルは、棋士レーティングを自己相関モデルで推定するものです。

完成したStanのモデルです。dataの構造はMemorandumを参考にしています。

model <- "
data {
  int N_kishi;
  int N_game;
  int N_year;
  int Winner[N_game];
  int Loser[N_game];
  int Year[N_game];
  int Career[N_kishi, 2];
}

parameters {
  real i_skill[N_kishi];
  real r_skill[N_kishi, N_year-1];
  real<lower=0> a;
  real<lower=0> s;
}

transformed parameters {
  real skill[N_kishi, N_year];
  real log_p[N_game];

  for(k in 1:N_kishi) {
    for (y in 1:Career[k,1])
      skill[k, y] = i_skill[k];

    for(y in (Career[k,1]+1):Career[k,2])
      skill[k, y] = skill[k, y-1] + r_skill[k, y-1];

    for (y in (Career[k,2]+1):N_year)
      skill[k, y] = skill[k,Career[k,2]];
  }

  for (g in 1:N_game)
    log_p[g] = log_inv_logit(a * (skill[Winner[g], Year[g]] - skill[Loser[g], Year[g]]));
}

model {
  for (k in 1:N_kishi) {
    target += normal_lpdf(i_skill[k] | 0, 100);
    for (y in Career[k,1]:(Career[k,2]-1))
      target += normal_lpdf(r_skill[k, y] | 0, s);
  }

  target += sum(log_p);

  target += uniform_lpdf(s | 0, 100);
}
"

モデル中のYearは公式戦準拠で年度単位となっています。スキル値も同様です。 全棋士の初期skillの事前分布は平均ゼロ、分散100の正規分布としました。 skillの自己相関モデルの標準偏差sの事前分布はゼロから100までの一様分布としました。 Careerは四段昇段以前の対局を含んだ、各棋士の実績の開始年度、終了年度のインデックスです。Career前のskillは初期skillと同じ値、Career後のskillはCareer終了年度のskillと同じ値にします。Career外の年度のskillは本来不要なので、サンプリング後に別処理でNAにしました。

サンプリングの実施

skill、a、sを収集対象パラメータとして、iter=12000、warmup=2000、thin=5、chains=4でサンプリングしました。対局結果データは収集した全データを使用しました。とてもものすごく時間がかかりました。よくやったもんです。

結果

Stanが計算してくれるRhatは全パラメータで1.1以下でした。収束したようです。やりました。

a, sの推定結果

skill差にかける係数aの平均は0.00671でした。イロレーティングの係数1/(400 * log(10))=0.0058と近いと言えるかもしれません。

meanse_meansd2.5%25%50%75%97.5%n_effRhat
a 0.00670692.7922e-050.00044769 0.0058647 0.0063853 0.0066965 0.0070188 0.0075937257.07 1.0092
s21.60418001.2963e-011.9210991818.245761520.184668821.536526222.904853925.5018467219.62 1.0108

skill差(レート差)と勝率の関係を可視化しました。イロレーティングとよく似ていますが偶然でしょう。適当に決めた初期スキルなどの事前分布が偶々こうなるような事前分布であった、ということだと思います。

skill差と勝率

skillの推定結果

skillの推定結果の統計量をこちらに載せています。

佐藤天彦のskill値の事後分布

Rhat的にはskill値はよく収束しているそうですが、安全安心のために事後分布を可視化してみます。実力制第十三代名人 佐藤天彦のskill値を年度ごとchainごとに可視化しました。

佐藤天彦

上下のヒゲの伸び方に差はありますが、それ以外はchain間の差がほとんど見られません。いい感じです。

棋士全年度についてskill値の平均値と中央値を比較すると大体同じでした。skill値の平均値と中央値の差の統計量を下に示します。どのskillも釣鐘型の綺麗な分布をしているんじゃないでしょうかと思います。

 diff_mean_median
 Min.   :-3.383  
 1st Qu.:-0.668  
 Median :-0.222  
 Mean   :-0.269  
 3rd Qu.: 0.172  
 Max.   : 1.734  

以降はskill値の代表値として平均値を使っていきます。

棋士レーティングとの比較

推定された2017年度のskill値(平均と95%区間)と棋士レーティング(2017年7月30日時点)を比較しました。Memorandumでも同様の可視化をしています。どうしても巨大な画像になりますね。

レーティングとの比較

大体整合的ですね。藤井聡太、大橋貴洸、大橋貴洸、杉本和陽の新人四人は95%区間が大きいです。それ以外の棋士の95%区間は大体同じような大きさに見えます。

下表は2017年度のskill値上位20人です。微差ながら藤井聡太佐藤天彦を抑えて堂々1位となりました。青嶋未来や近藤誠也もskill値では評価高く、まだレーティングが追いついていない状態と言えます。

skill順位名前skill平均skill95%区間skill95%区間レーティングレーティング順位
1 藤井聡太 318.15 206.71 438.42 1703 32
2 佐藤天彦 317.40 243.53 397.99 1884 1
3 羽生善治 315.03 243.92 394.22 1849 4
4 豊島将之 313.00 240.05 390.72 1878 2
5 渡辺明 288.24 213.38 367.98 1828 7
6 菅井竜也 286.98 216.12 365.51 1856 3
7 稲葉陽 274.81 201.85 354.58 1837 6
8 久保利明 274.77 202.91 352.13 1839 5
9 永瀬拓矢 258.59 186.67 334.47 1825 8
10 糸谷哲郎 252.53 182.40 327.20 1808 9
11 千田翔太 244.15 173.17 319.49 1790 11
12 斎藤慎太郎240.46 169.81 313.21 1801 10
13 佐々木勇気238.27 170.34 310.52 1778 13
14 広瀬章人 234.33 163.60 309.05 1780 12
15 山崎隆之 229.63 160.52 304.31 1776 14
16 中村太地 226.65 153.89 300.82 1775 15
17 三浦弘行 224.80 151.06 303.74 1767 17
18 青嶋未来 224.12 147.35 302.01 1715 26
19 近藤誠也 224.00 142.70 307.73 1701 33
20 松尾歩 219.55 148.18 295.29 1766 18

棋士のskill値の推移

せっかく50年余りの推定を行ったので、全棋士のskill値の平均の推移をプロットしました。

15MBのpngファイル

デカイです。棋士名が重なって見づらいです。個別の棋士の推移を追うことができないですね。

全体の傾向として、どうもskill値の分布の裾が時間経過とともに広がっていることがわかります。 試しに5年単位でskill値の平均のバイオリンプロットを描いてみました。

skillの密度分布の推移

分布の中心とばらつきが一定だと年度間でのskill値を容易に比較できて大変便利だと思う、のですが。

  1. 中央値、平均値とも上昇傾向, skill値はインフレしている
  2. ばらつきも大きくなる傾向があるような気がするが、はっきりしない。1995年付近を境界にしてばらつきが大きくなった=棋力の格差が広がった、ように見える。

これでは、年度をまたいだskillの値そのものの比較はできませんね。具体的に今回のモデルの何が悪いのかは分かりません。

ちなみにイロレーティングもインフレ・デフレするそうです。

歴代のトップ棋士

各年度におけるskill値の順位TOP10を可視化しました。黒丸は公式戦デビュー時点でランクインしたことを示します。黒字のラベルは翌年度にランク外に陥落したことを示します。

トップ棋士のskill推移

大山→中原→羽生と続く棋界の頂点の推移がよくわかりますね。また、中原誠羽生善治とも、デビュー時で2位にランクインし、すぐに1位に上り詰めています。先程藤井聡太がskill1位すごいすごいと書いたわけですが、確かにすごいが前例はあるのでした。羽生世代やポスト羽生世代には、デビュー即トップ10入りという棋士が結構います。もちろんデビュー直後なので信頼区間の幅も大きいわけですが。

谷川浩司渡辺明は一度も1位になっていません。これは違和感があります。谷川浩司については、中原以降羽生以前に覇権を握ったイメージがあったのですが、今回のモデルでは、skill上はそのような時代は無かった、となります。

羽生世代が80年代後半にドカドカと現れて、2010年台にドカドカと退場する、世代交代の様が一目瞭然です。羽生善治も2016年度に1位の座を明け渡したようです。レーティングもそんな感じですね。次なる頂点の候補は佐藤天彦藤井聡太と、最近の若手では珍しくデビュー時点でトップ10入りの菅井竜也でしょうか。私は糸谷哲郎推しですが、この人ポカが多いんですよね。

山田道美、村山聖はトップ10のまま夭折しました(村山聖は最期は休場)。山田道美は山田定跡の考案者ですね。大山康晴相手にタイトルもとりました。棋士番号制定前の棋士であるためか、将棋連盟棋士データベースには個別ページさえもありません。

TOP10の棋士のスキルはおおむね右肩上がりに上昇しています。その中で羽生善治は、1995年から2010年まで、2位の棋士相手に常に50以上のskill差をつけています。やはり化物ですね。

棋士の全盛期と年齢

棋士の全盛期と年齢の関係を可視化しました。年齢はデータ収集のところで書いた生年月日をもとに、年度ベースの数え歳で計算しました。つまり平成27年度生まれは平成27年度に1歳、という換算になります。

全盛期をskill値そのもので求めるのは危ないので、期歴におけるskill値の年度順位が最高の年度を全盛期としました。例えば羽生善治は15歳から46歳まで1位を32年間維持したので、この間ずっと全盛期という扱いになります。また、この計算で全盛期が年度期間の両端(1964年、1965年、2016年、2017年)となった棋士は、真の全盛期は計算期間の範囲外にあるものと思われるため、除外しました。

棋士の全盛期年齢

棋士全体ではskill値は20代半ばに全盛期が来ることが多いようです。トップ棋士の場合は山がもっとのっぺりしてます。

40過ぎで全盛期が来て(あるいは全盛期を維持して)、それがトップ20以内の棋士は以下になります。花村元司以上は、1963年以前に全盛期があった可能性があると思います。加藤一二三は名人獲得直前が全盛期ですね。加藤一二三は早熟の天才というイメージがありましたが、今回収集したデータの範囲では40歳前後で順位を上げており、実に不思議な棋士です。「大山康晴に苛められすぎたのが中年になって吹っ切れた」と、中原誠渡辺明との対談で言っていたような記憶があります(たしか将棋世界誌上)。

idname年度skill値の平均順位誕生年度年齢
1026 大山 康晴 1966 197.1745 1 1922 45
1035 原田 泰夫 1966 -4.2594 19 1922 45
1035 原田 泰夫 1967 -5.3930 19 1922 46
1039 花村 元司 1966 7.0797 16 1917 50
1064 加藤 一二三1981 131.3172 3 1939 43
1085 米長 邦雄 1982 141.3653 2 1943 40
1175 羽生 善治 2009 337.4900 1 1970 40
1175 羽生 善治 2010 349.6142 1 1970 41
1175 羽生 善治 2011 353.2264 1 1970 42
1175 羽生 善治 2012 358.3190 1 1970 43
1175 羽生 善治 2013 352.0057 1 1970 44
1175 羽生 善治 2014 344.3521 1 1970 45
1175 羽生 善治 2015 330.5932 1 1970 46
1204 三浦 弘行 2015 230.0110 13 1973 43

今回のモデルは単純なので、推定されたskillの平均の推移は、1963年からイロレーティングを計算した場合と大きな違いはないと思います。 デビュー直後の棋士の強さを推定できるのことが優位点でしょうか。 モデルにもっと棋士個人差のパラメタを組み込むと面白いと思います。例えばskillの自己相関の標準偏差や、勝負ムラなど。

色々な可視化ができたのは楽しかったです。

おわり。


  1. 久保拓弥著 「データ解析のための統計モデリング入門」という本の読書会資料をググって漁っていたところ、知りました。

やってみた。ドラゴンボールの戦闘力推定

ドラゴンボールの勝敗結果から戦闘力を推定する, StatModeling Memorandumという実験記事があります。面白そうなのでやってみました。

準備

RとStanを使いました。tidyrとdplyrはRで使われる変な名前のデータ整形ライブラリです。

set.seed(1)

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(tidyr)
library(dplyr)

データ

games.txtは戦闘データ、senshis.txtはキャラクタ名と戦闘力の文献値が入っています。値はMemorandumの記事から拝借しました。

%>%はtidyrやdplyrで使われる、メソッドチェイン的な構文が使える演算子です。joinを使って戦績にsenshi.txtを連結した結果を出しています。

games <- read.table("games.txt", header=T)
senshis <- read.table("senshis.txt", header=T)

games %>%
    left_join(senshis, c("winner"="ID")) %>%
    left_join(senshis, c("loser"="ID"), suffix = c(".winner", ".loser"))
winnerlosername.winnerpower.winnername.loserpower.loser
1 2 ベジータ18000 悟空 8000
1 3 ベジータ18000 ナッパ 4000
1 5 ベジータ18000 栽培マン1200
2 3 悟空 8000 ナッパ 4000
3 4 ナッパ 4000 クリリン1770
4 5 クリリン 1770 栽培マン1200

モデル

モデルです。元記事のBUGSのモデルをまったくそのままStanに書き直したつもりです。

model <- "
data {
  int N_member;
  int N_game;
  int Winner[N_game];
  int Loser[N_game];
}
parameters {
  real winner_p[N_game];
  real loser_p[N_game];
  real power[N_member];
}
model {
  for (game in 1:N_game) {
    target += normal_lpdf(winner_p[game] | power[Winner[game]], 1);
    target += normal_lpdf(loser_p[game] | power[Loser[game]], 1);
    target += bernoulli_lpmf(1 | step(winner_p[game] - loser_p[game]));
  }
  for (member in 1:N_member)
    target += normal_lpdf(power[member] | 0, 100);
}
"

実行と結果

stan関数でサンプリングを実行します。初期値を与えないとサンプリングに失敗しました。とはいえ下で使用したinit関数は適当すぎるかもしれませんが、とりあえず気にしないことにします。

N_member <-  nrow(senshis)
N_game <- nrow(games)

data <- list(
  N_member = N_member,
  N_game   = N_game,
  Winner   = games$winner,
  Loser    = games$loser
)

init <- function(...) {
  list(
    winner_p = rep(1, N_game),
    loser_p  = rep(0, N_game),
    power    = rep(0, N_member)
  )
}

chains <- 4
 
fit <- stan(
  model_code = model,
  data       = data,
  init       = lapply(1:chains, init),
  pars       = c('power'),
  iter       = 22000,
  warmup     = 2000,
  thin       = 10,
  chains     = chains
)

結果です。powerしか見ていませんが、Rhatが1.1以下なのでよく収束していると言えそうです。

print(fit)
Inference for Stan model: 174a3ef46412238a4997765f45f7964d.
4 chains, each with iter=22000; warmup=2000; thin=10; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

            mean se_mean    sd    2.5%     25%     50%    75%  97.5% n_eff Rhat
power[1]  117.90    3.06 65.08    0.42   73.67  113.45 157.94 258.26   452 1.00
power[2]   52.57    2.89 55.15  -53.83   14.25   52.60  89.81 165.15   365 1.01
power[3]    5.47    2.82 54.18 -100.56  -32.49    7.71  42.99 106.10   370 1.01
power[4]  -44.23    2.64 55.49 -156.75  -82.17  -42.78  -5.11  57.88   443 1.00
power[5] -110.90    2.85 65.68 -251.55 -152.42 -105.90 -63.64   5.47   531 1.00
lp__      -47.05    0.07  2.87  -53.58  -48.74  -46.74 -45.00 -42.42  1850 1.00

Samples were drawn using NUTS(diag_e) at Sun Jul 30 23:08:37 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

ただしサンプリング時に以下の警告が出ます。意味はわかりませんがモデルが良くないと言われているような気がします。試しにadapt_deltaを増やしたサンプリングもやってみましたが、エラーは消えませんでした。これも今回は気にしないことにします。

 警告メッセージ:
1: There were 5873 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
2: There were 520 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
3: Examine the pairs() plot to diagnose sampling problems

可視化

さて楽しい楽しい可視化の時間です。むしろ可視化がしたいがためにこの記事を書いたようなものです。

ggplot2というRではメジャーな可視化ライブラリを使って、Memorandumの確率密度プロットを再現しました。バイオリンプロットと呼ぶそうです。

まずstanのサンプリング結果オブジェクトからサンプリング値にアクセスします。extract関数を使いますが、extractはtidyrロード時に同名の別関数でマスクされているので、::でnamespaceを指定して呼び出します。strは文字列変換ではなく、オブジェクトの構造を出力する関数です。stringではなくstructureです。

str(rstan::extract(fit, pars="power"))
List of 1
 $ power: num [1:8000, 1:5] 115.5 141.9 125.3 104.5 57.7 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ iterations: NULL
  .. ..$           : NULL

可視化用にデータを加工します。戦闘力の文献値は常用対数にします。

senshis$power <- log10(senshis$power)

D <- rstan::extract(fit, pars="power")$power %>%
  data.frame %>%
  setNames(senshis$name) %>%
  gather(name, power) %>%
  inner_join(senshis, c("name"="name"), suffix=c('.est.','.lit.'))

D$name = factor(D$name, levels=senshis$name)
head(D)

CI <- D %>%
  group_by(name) %>%
  summarise(power.lit.=mean(power.lit.),
            median=median(power.est.),
            lower_95=sort(power.est.)[length(power.est.)*0.025],
            upper_95=sort(power.est.)[length(power.est.)*0.975])
CI
namepower.est.IDpower.lit.
ベジータ 115.497261 4.255273
ベジータ 141.919341 4.255273
ベジータ 125.329281 4.255273
ベジータ 104.522831 4.255273
ベジータ 57.720201 4.255273
ベジータ 53.393671 4.255273
namepower.lit.medianlower_95upper_95
ベジータ 4.255273 113.446818 0.04747007258.253949
悟空 3.903090 52.597184 -53.84137102165.144239
ナッパ 3.602060 7.713342 -100.58709283106.098659
クリリン 3.247973 -42.779770 -157.13586821 57.882540
栽培マン 3.079181 -105.896344 -251.56942070 5.471327

ggplot2で可視化します。まったく慣れていないので大変でしたが、ググれば情報は出て来るのでやりようはあります。

p <- ggplot(D, aes(power.lit., power.est., fill=name)) +
  geom_violin(size=0, adjust=1.5) +
  geom_pointrange(CI, mapping=aes(x=power.lit., y=median, ymin=lower_95, ymax=upper_95, colour=name), size=1.5) +
  scale_fill_manual(values=rep("grey70", 5)) +
  scale_y_continuous(limits=c(-500, 500), breaks=seq(-500, 500, 250)) +
  xlab("log10(Literature Power Level)") +
  ylab("Estimated Power Level") +
  theme(plot.title=element_text(size=18)) +
  theme(axis.text.x=element_text(size=14)) +
  theme(axis.text.y=element_text(size=14)) +
  theme(axis.title.x=element_text(size=18)) +
  theme(axis.title.y=element_text(size=18)) +
  theme(legend.title=element_text(colour="black",size=18)) +
  theme(legend.text=element_text(colour="black",size=18))

png(filename='result.png', width=700, height=600)
plot(p)
garbage <- dev.off()

可視化

Memorundomと同じようなグラフがかけました。ベジータの上方向と栽培マンの下方向が伸び伸び、ベジータ-悟空間、ナッパ-クリリン間の戦闘力差が小さいなど、特徴もよく似ています。Stanの結果も可視化処理も問題ないようです。

おわり。次は階層ベイズモデルで勝敗データからプロ棋士の強さを推定する, StatModeling Memorandumをやってみたいです。

ギブスサンプリングでベイズ推定

ギブスサンプリングは確率分布からのサンプリング手法の一つ、各変数について条件付き確率さえわかればサンプリングできる。ベイズ推定での事後分布サンプリングに使われることが有る(Just another Gibbs samplerというそのものズバリなベイズ推定のオープンソースプロダクトもある)。

以下はギブスサンプリングに関する参考文献。2は可視化もあって大変分かりやすい。直感的に理解できるのがギブスサンプリングのいいところだと思う。

  1. Gibbs_sampling, wikipedia
  2. 可視化で理解するマルコフ連鎖モンテカルロ法(MCMC), ほくそ笑む

ここではギブスサンプリングを実装してベイズ推定を実施する。ギブスサンプリングは簡単なアルゴリズムなので、多分、実装も簡単であろうと思う。

実装は以下(Campbell)を参考にした。これを読めば誰でもギブスサンプリングでベイズ推定ができる、とても分かりやすい記事。

Gibbs sampling for Bayesian linear regression in Python, Kieran Campbell

テストデータ、モデルは以下の記事(Memorandum)のものをそのまま利用した。記事中のRstanのコードは推定結果の答え合わせに利用した。このブログは好奇心をそそられる記事に溢れており、眺めるだけで賢くなった気分になれる。

WAICとWBICを事後分布から計算する, StatModeling Memorandum

データと統計モデル

2成分混合正規分布から100点を生成しテストデータとする。MemorandumがRを使っているので、それにならう。

from subprocess import Popen, PIPE

R = Popen(["R", "--vanilla", "--silent"], stdin=PIPE, stdout=PIPE)

R.stdin.write(b'''
set.seed(1)
N      <- 100
a_true <- 0.4
mean0  <- 0
mean1  <- 3
sd0    <- 1
sd1    <- 1
Y      <- c(rnorm((1-a_true)*N, mean0, sd0), rnorm(a_true*N, mean1, sd1))

write.table(Y, file="points.csv", sep=",", row.names=F, col.names=F)
q()
''')
R.communicate()
R.kill()

統計モデルは2つ。それぞれについてベイズ推定を行う。

  1. 単一の正規分布モデル
  2. 2成分混合正規分布モデル

Pythonセットアップ

CampbellにならってPython3を用いる。pandasというのはRのdata.frame的なものらしい。多彩なデータ加工ができるようだが趣味人には使い方が分からない。Anacondaは怖いので使わない。

$ python -m venv .venv
$ source .venv/bin/activate
$ pip install numpy scipy pandas seaborn

以下はpythonの共通設定的なもの。

import numpy as np
np.random.seed(0)

import scipy as sp

import pandas as pd
pd.set_option('display.width', 200)

import seaborn as sns
from seaborn import plt
plt.rcParams['figure.figsize'] = (12, 6)

可視化のテスト、R経由で作成したテストデータを可視化してみる。

Y = pd.read_csv("points.csv", header=None)
p = sns.distplot(Y, bins=20)
plt.savefig("points.png")

f:id:kazufusa1484:20170621132727p:plain

混合分布ってのはこういう形をしているようだ。

ギブスサンプリングでベイズ推定のおさらい

Campbellの記事でギブスサンプリングを復習。

今手元に N 個のモデルパラメタ  \theta_j, j=1,\ldots,n とデータ  y があり、ここから事後分布  p(\theta_j \mid y) を得たいとする。ギブスサンプリングで事後分布をサンプリングするには、条件付き確率 p(\theta_j\mid\theta_1, \ldots, \theta_{j-1}, \theta_{j+1}, \ldots, \theta_n, y) を用いて以下の手順をとる。

  1. 適当な初期値 \theta_j^{(i)} を与える.
  2. すべての jについて  \theta_j^{(i+1)} \sim p(\theta_j\mid\theta_1^{(i+1)},\ldots, \theta_{j-1}^{(i+1)},\theta_{j+1}^{(i)}, \ldots, \theta_n^{(i)}, y) をサンプリング.
  3. 2.を予め決めたイテレーション数だけ繰り返す.

イテレーション数が十分であれば、各パラメタのサンプリング結果の密度分布を、各パラメタの周辺化した事後分布として扱うことができる、らしい。

確率モデルの実装(モデル1: 単一の正規分布モデル)

モデルの概要

ここではデータ  y が従う正規分布の平均  \mu 標準偏差  \sigma を推定する。


y_i \sim \mathcal{N}(\mu, \sigma^2), i=1,\ldots, N

データ  y の条件付き確率。


p(y_1, \ldots, y_N \mid \mu, \sigma) = \prod_{i = 1}^{N} \mathcal{N}(y_i\mid \mu, \sigma^2)

パラメタには事前分布は、簡単のため無情報事前分布とする。


\mu \sim \rm{uniform}(-\infty, \infty)


\sigma \sim \rm{uniform}(0, \infty)

このときパラメタの確率は一定になる。


 p(\mu)=\rm {constant}


 p(\sigma) = \rm {constant}

 \mu のサンプリング

何はなくとも条件付き確率がわからなければサンプリングできない。 定数項と \mu に独立な項を消去する。



\begin{eqnarray}
p(\mu \mid \sigma, y) &=& \frac {p(\mu, \sigma, y)} {p(\sigma, y)} \\
&\propto& p(y \mid \mu, \sigma) p(\mu) p(\sigma) \\
&\propto& p(y \mid \mu, \sigma) \\
\end{eqnarray}

データ  y の条件付き確率が出てきた。



\begin{eqnarray}
p(y \mid \mu, \sigma) &=& \prod_{i = 1}^{N} \mathcal{N}(y_i\mid \mu, \sigma^2) \\
 &=&  \prod_{i = 1}^{N} \left(\frac {1} {\sqrt{2\pi\sigma^2}} exp \left(- \frac {(y_i- \mu)^2} {2\sigma^2} \right) \right)\\
 &\propto&  \prod_{i = 1}^{N} exp\left(-\frac {\mu^2} {2\sigma^2} + \frac {y_{i} \mu} {\sigma^2}\right)\\
 &=& exp\left(-\frac {N} {2\sigma^2} \mu^2 + \frac { \sum_{i = 1}^{N} y_{i}} {\sigma^2} \mu \right)\\
 &\propto& \frac {1} {\sqrt{2\pi \frac{\sigma^2} {N}}} exp \left(- \frac { \left(\mu - \frac {1} {N} \sum_{i = 1}^{N} y_{i} \right)^2} {2 \frac{\sigma^2} {N}} \right)        \\
 &=&  \mathcal{N}\left( \frac{1} {N} \sum_{i=1}^{N} y_i, \frac {\sigma^2} {N} \right) \\
\end{eqnarray}

最後の二行で、規格化されていない尤度関数からうまいこと正規分布が得られた。

これはつまり、ある確率変数  x の尤度が exp(x の二次式) であるとき、確率変数  x 正規分布をとる、ということである。対数尤度が x の二次式なら正規分布!、と憶えるとよさそう。

このようにして  \mu の条件付き確率がわかった。


 \mu \mid \sigma, y \sim \mathcal{N}\left( \frac{1} {N} \sum_{i=1}^{N} y_i, \frac {\sigma^2} {N} \right)

 \mu のサンプリングコード、正規分布からのサンプリングなのでとても簡単。

def sample_mu(y, N, s):
    mean = np.sum(y) / N
    variance = s * s / N
    return np.random.normal(mean, np.sqrt(variance))

 \sigma のサンプリング

つづいて  \sigma  \mu と同様に条件付き確率の比例成分として  y の事後分布が得られる。このままでは計算しづらいので、対数尤度  lp(\sigma) を解き、 \sigma がどのような分布に従うか調べる。


 p(\sigma \mid \mu, y) \propto p(y \mid \mu, \sigma)


\begin{eqnarray}
lp(\sigma) &=& Nlog\left( \frac {1} {\sqrt{2\pi\sigma^2}}\right) -  \frac {1} {2\sigma^2} \sum_{i=1}^N \left(y_i- \mu \right)^2 \\
&=& -Nlog\sigma - \frac {1} {2\sigma^2} \sum_{i=1}^N \left(y_i- \mu \right)^2 - \frac {N} {2} log 2\pi\\
\end{eqnarray}

対数尤度が \sigma の二次式ではないので、とりあえず正規分布ではない。ここでガンマ分布(Gamma distribution:wikipedia)を導入する。

ガンマ分布は  \alpha\beta の2パラメタで定義される分布であり、確率変数  x がガンマ分布に従うとき確率密度は  p(x\mid \alpha, \beta) \propto \beta^\alpha x^{\alpha-1}e^{-\beta x} となる。 x の依存項のみを抜き出した対数尤度 lp(x) は以下になる。


 lp(x\mid \alpha, \beta) = (\alpha - 1)logx - \beta x

 lp(\sigma) と似たような形をしている。ここで、 \sigma から \tau = 1/\sigma^{2} \tau は精度と呼ばれるもの)への変数変換を行い、


 lp(\sigma) = \frac{N} {2} log\tau -  \frac {\tau} {2} \sum_{i=1}^N \left(y_i- \mu \right)^2- \frac {N} {2} log 2\pi

より、 \tau の条件付き確率は以下のガンマ分布に従うことが示せた。


 \tau \mid \mu, y \sim \rm {Gamma} \left(\frac{N} {2}+1, \frac {\sum_{i=1}^N \left(y_i- \mu \right)^2} {2}  \right)

サンプリングコード、 \tau をサンプリングし、後に  \sigma に変換する、という実装。

def sample_s(y, N, mu):
    alpha = N / 2 + 1
    residuals = y - mu
    beta = np.sum(residuals * residuals) / 2
    tau = np.random.gamma(alpha, 1 / beta)
    return 1 / np.sqrt(tau)

サンプリングの実施

条件付き確率が得られたので、Campbellを参考にサンプラを実装する。

イテレーションごとに 各  y の対数尤度(log_likelihoods)を計算している。これはMemorandumを参考にしたもので、WAICというモデルの汎化性能を評価する指標があり、その算出に必要。ここではMemorandumの結果と一致することを確認するために計算している。

def model1(y, iters, init):
    mu = init["mu"]
    s = init["s"]
    N = len(y)
    
    trace = np.zeros((iters, 2))
    log_likelihoods = np.zeros((iters, N))
    
    for i in range(iters):
        mu = sample_mu(y, N, s)
        s = sample_s(y, N, mu)
        trace[i, :] = np.array((mu, s))
        
        norm = sp.stats.norm(mu, s)
        log_likelihoods[i, :] = np.array([np.log(norm.pdf(x)) for x in y])
    
    trace = pd.DataFrame(trace)
    trace.columns = ['mu', 's']
    
    log_likelihoods = pd.DataFrame(log_likelihoods)
    
    return trace, log_likelihoods

初期値。

init = {"mu": np.random.uniform(-100, 100), "s": np.random.uniform(0, 100)}
print(init)
{'mu': 9.762700785464943, 's': 71.51893663724195}

サンプリング実施。

iters = 5000
y = np.array([float(x) for x in open("points.csv", "r").read().strip().split("\n")])
trace, log_likelihoods = model1(y, iters, init)

結果

サンプリング結果をイテレーションを横軸としてプロットしたものをトレースプロットと呼ぶ。これが一定の範囲内で上下ふらふらしていればひとまずサンプリングは成功と言える。見た感じは良さそう。1

traceplot = trace.plot()
traceplot.set_xlabel("Iteration")
traceplot.set_ylabel("Parameter value")
traceplot.get_figure().savefig("model1_trace.png")

f:id:kazufusa1484:20170621132822p:plain

パラメタごとのサンプリング結果のヒストグラム。サンプリング結果の後半だけをtrace_burntに切り出している。ギブスサンプリングを含めたマルコフ連鎖モンテカルロサンプリングでは、サンプリングの序盤は暖機運転扱いで、結果から除外するのが通例であるらしい2

trace_burnt = trace[int(len(trace)/2):]
hist_plot = trace_burnt.hist(bins = 30, layout = (1,2))
traceplot.get_figure().savefig("model1_hist.png")

f:id:kazufusa1484:20170621132840p:plain

MemorandumのRコードからWAIC算出関数を移植。

def waic(log_likelihoods):
    training_error = -np.log(np.exp(log_likelihoods).mean(axis=0)).mean()
    functional_variance_div_N = (np.power(log_likelihoods, 2).mean(axis=0) -
                                 np.power(log_likelihoods.mean(axis=0),2)).mean()
    return training_error + functional_variance_div_N

統計量等と合わせて表示する。WAICはMemorandumでは1.980なので、よく一致している。成功したらしい。

print(trace_burnt.describe().T)
print("waic:", waic(log_likelihoods))
print("np.mean(y):", np.mean(y))
print("np.std(y):", np.std(y))
     count      mean       std       min       25%       50%       75%       max
mu  2500.0  1.307089  0.173329  0.596588  1.193417  1.307739  1.419458  1.945270
s   2500.0  1.728932  0.119035  1.411500  1.644560  1.725093  1.802547  2.212252
waic: 1.98051651897
np.mean(y): 1.30888736691
np.std(y): 1.7214151253

確率モデルの実装(モデル2: 2成分混合モデル)

次のモデルに進む。2つの正規分布の混合モデル。

モデルの概要

何も考えずにモデルを書くと下になる、と思う。


y_i \sim (1-a)  \mathcal{N}(\mu_0, \sigma_0) + a \mathcal{N}(\mu_1, \sigma_1)

ところがこのモデルはどういじっても条件付き確率が導出できない。どうも正規分布の足し算が難易度高いらしい。無理。

困ったのでmixture、gaussian、gibbs等で適当にググると、混合分布モデルの場合はデータ点が属するカテゴリを推定するのがよい、と書いてあるらしい難しい内容の記事・文献がいくつか出てた3。要するにデータ  y_i の属する正規分布カテゴリを表す確率変数  z_i を導入し、


 z_i \sim \text {Bernoulli} (a)

y_i \mid z_i \sim  \mathcal{N}(\mu_{z_i}, \sigma_{z_i}^{2})

とするといいらしい。Bernoulliは確率aで1、(1-a)で0となる離散確率分布、曲がったコインのトスみたいなもの。これはつまり、今回のテストデータ作成をそのままモデル化したのと同じことである。 y の条件付き確率は、


p(y_1, \ldots, y_N \mid \mu_1, z_1, \ldots, z_N) = \prod_{i = 1}^{N} \mathcal{N}(y_i\mid \mu_{z_i}, \sigma_{z_i}^{2})

となり、掛け算なので取り扱いが楽そう。確かに先に進めそうである。

 \mu_0 = 0, \sigma_0 = 1, \sigma_1 = 1 は既知とする。 \mu_1  a の事前分布は以下とする。


\mu_1 \sim \text{uniform}(-\infty, \infty)


a \sim \text{uniform}(0, 1)

 \mu_1 のサンプリング

 \mu_1 の条件付き確率。



\begin{eqnarray}
p(\mu_1 \mid a, z, y) &\propto& p(y \mid \mu_1, z) \\
&=& \prod_{i = 1}^{N} \mathcal{N}(y_i \mid \mu_{z_i}, \sigma_{z_i}^{2}) \\
&\propto& \prod \mathcal{N}(y_i\mid \mu_{z_i}, \sigma_{z_i}^{2})\mid_{z_i=1} \\
\end{eqnarray}

考え方はモデル1の  \mu と同じで、正規分布になる。ただしデータは  y 全体ではなく、 \mathcal{N}( \mu_1, 1) 正規分布に属するデータ( y_i \mid_{z_i=1} )のみを用いる。


 \mu_1 \mid z_1, \ldots, z_N, y \sim \mathcal{N}\left( \frac{1} {\sum^{N}_{i=1} z_i} \sum y_i\mid_{z_i=1}, \frac {1} {\sum^{N}_{i=1} z_i} \right)

サンプリングコード。コード中のzは  z_i の配列だが、要素が0か1なので、  y_i\mid_{z_i=1} の総和や総数は以下のように楽ちんに計算できる。

def sample_mu_1(y, z):
    N1 = np.sum(z)
    mean = np.sum(y * z) / N1
    variance = 1 / N1
    return np.random.normal(mean, np.sqrt(variance))

 a のサンプリング

a の条件付き確率。 z_i は0か1の離散値なので、総和を取ると \mathcal{N} (\mu_1, 1) に属するデータ点の総数になる。



\begin{eqnarray}
p(a \mid \mu_1, z, y) &\propto& \prod_{i=1}^{N} \text{Bernoulli}(z_i \mid a) \\
&=& (1-a)^{N-\sum_{i=1}^{N} z_i} a^{\sum_{i=1}^{N} z_i}
\end{eqnarray}

例によって上の二行目は規格化すると名前のある分布になるのだろうと期待される。ベータ分布  \left( p(x) = \frac {x^{\alpha-1}(1-x)^{\beta-1}} {B(\alpha, \beta)}\right) であるようだ。



p(a\mid \mu_1, z, y) = \text {Beta} \left( \sum_{i=1}^{N} z_i + 1, N - \sum_{i=1}^{N} z_i + 1 \right)
def sample_a(z, a):
    alpha = np.sum(z) + 1
    beta = len(z) - alpha + 2
    return np.random.beta(alpha, beta)

 z_i のサンプリング

 z_i の条件付き確率。 z_i は0 or 1 の離散値なので、尤度さえ分かればサンプリング可能。



\begin{eqnarray}
p(z_i \mid \mu_1, a, z_1, \ldots,z_{i-1}, z_{i+1}, \ldots, z_N, y) &\propto& p(y \mid \mu_1, z_1, \ldots, z_N) p(z_i\mid a) \\
&\propto& \mathcal{N}(y_i \mid \mu_{z_i}, \sigma_{z_i}^{2}) \times \text{Bernoulli}(z_i \mid a) \\
\end{eqnarray}

コードでは一つの関数で全  z_i をサンプリングしている。 \mathcal{N}(0, 1) は使いまわせるのでグローバルに定義している。 zのサンプリング結果を返すのではなく、zの配列を更新しているため、関数名がこれまでと異なる。

norm0 = sp.stats.norm(0, 1)

def sample_and_update_z(y, mu_1, a, z):
    norm1 = sp.stats.norm(mu_1, 1)

    for i in range(len(y)):
        d0 = (1-a) * norm0.pdf(y[i])
        d1 = a * norm1.pdf(y[i])
        z[i] = 0 if np.random.uniform() * (d0+d1) < d0 else 1

サンプリングの実施

サンプリング実行コードと初期値。

def model2(y, iters, init):
    mu_1 = init["mu_1"]
    a = init["a"]
    z = init["z"]
    N = len(y)
    
    trace = np.zeros((iters, 2))
    log_likelihoods = np.zeros((iters, N))
    
    for i in range(iters):
        mu_1 = sample_mu_1(y, z)
        a = sample_a(z, a)
        sample_and_update_z(y, mu_1, a, z)
        trace[i, :] = np.array((mu_1, a))
        
        norm1 = sp.stats.norm(mu_1, 1)
        log_likelihoods[i, :] = \
            np.array([np.log((1-a)*norm0.pdf(x) + a*norm1.pdf(x)) for x in y])
    
    trace = pd.DataFrame(trace)
    trace.columns = ['mu_1', 'a']
    
    log_likelihoods = pd.DataFrame(log_likelihoods)
    
    return trace, log_likelihoods
init = {"mu_1": np.random.uniform(-100, 100), "a": np.random.rand(), "z": np.random.randint(0, 2, 100)}
print(init)
{'mu_1': -99.69811311464245, 'a': 0.7177148886002629, 'z': array([1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1,
       1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0,
       1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0,
       1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1,
       0, 1, 0, 1, 0, 0, 1, 1])}

サンプリングを実施。

iters = 5000
y = np.array([float(x) for x in open("points.csv", "r").read().strip().split("\n")])
trace, log_likelihoods = model2(y, iters, init)

結果

traceplot = trace.plot()
traceplot.set_xlabel("Iteration")
traceplot.set_ylabel("Parameter value")
traceplot.get_figure().savefig("model2_trace.png")

f:id:kazufusa1484:20170621132905p:plain

trace_burnt = trace[int(len(trace)/2):]
hist_plot = trace_burnt.hist(bins = 30, layout = (1,2))
traceplot.get_figure().savefig("model2_hist.png")

f:id:kazufusa1484:20170621132919p:plain

統計量とWAIC。Memorandumでは1.913とのこと。ほぼ同じである。

print(trace_burnt.describe().T)
print("waic:", waic(log_likelihoods))
       count      mean       std       min       25%       50%       75%       max
mu_1  2500.0  3.083692  0.205765  2.199600  2.940599  3.085343  3.223398  3.768049
a     2500.0  0.397030  0.055852  0.214459  0.360779  0.395766  0.433929  0.601889
waic: 1.91467884523

Rstanでの結果

推定結果は良さそうだし、WAICもMemorandumとほぼ同じだ。ついでにMemorandumに書かれたRstanコードを実行し、推定値の統計量を比較する。

from subprocess import check_output

print(check_output(["Rscript", "model2.r"]).decode("utf8"))
# results of the mixture normal distribution model with Rstan
        mean      se_mean         sd      2.5%      25%       50%       75%    97.5%    n_eff     Rhat
mu 3.0914330 0.0012293010 0.20128997 2.6958750 2.955486 3.0909330 3.2286982 3.485634 26811.91 1.000001
a  0.3965491 0.0003816836 0.05627476 0.2907154 0.357534 0.3955335 0.4343038 0.509682 21738.04 1.000274
WAIC: 1.914037

これもよく一致した。やはりうまくいったようだ。

上で使用したRコードはこちら。

options(width=200)
Y <- read.table("points.csv")
data   <- list(N=length(Y), Y=Y)

model2 <- "
data {
  int<lower=1> N;
  vector[N] Y;
}

parameters {
  real<lower=0, upper=1> a;
  real<lower=-50, upper=50> mu;
}

model {
  for(n in 1:N){
    target += log_sum_exp(
      log(1-a) + normal_lpdf(Y[n] | 0, 1),
      log(a) + normal_lpdf(Y[n] | mu, 1)
    );
  }
}

generated quantities {
  vector[N] log_likelihood;
  int index;
  real y_pred;
  for(n in 1:N)
    log_likelihood[n] = log_sum_exp(
      log(1-a) + normal_lpdf(Y[n] | 0, 1),
      log(a) + normal_lpdf(Y[n] | mu, 1)
    );
  index = bernoulli_rng(a);
  y_pred = normal_rng(index == 1 ? mu: 0, 1);
}
"

sink(file="/dev/null")
suppressMessages({
  library(rstan)
  fit <- stan(model_code=model2, data=data, iter=11000, warmup=1000, seed=123)
})
sink()
cat("# results of the mixture normal distribution model with Rstan\n")
print(summary(fit)$summary[c("mu", "a"), ])

waic <- function(log_likelihood) {
  training_error <- - mean(log(colMeans(exp(log_likelihood))))
  functional_variance_div_N <- mean(colMeans(log_likelihood^2) - colMeans(log_likelihood)^2)
  waic <- training_error + functional_variance_div_N
  return(waic)
}

cat(sprintf("WAIC: %f", waic(extract(fit)$log_likelihood)))

おわり。


  1. 見た感じでは許されない状況にある場合は、2つ以上のギブスサンプリングを回して、どれだけ結果が似ているかを評価する。Rhat(potential scale reduction factor)という指標がメジャー。簡単な統計量なのに実装はいくつかあってどれを使えばいいのかよく分からない。どれでもいいのかもしれない。

  2. 他に、サンプリングとサンプリングの間に何回かのイテレーションを挟むことがある。 目的は不明。たくさんイテレーションを回したいが、サンプリング数が多すぎると結果の処理が辛いため?

  3. 日本語だとベイズ混合モデルにおける近似推論③ ~崩壊型ギブスサンプリング~, 作って遊ぶ機械学習。 があった。内容は難しくて3%も理解できない。