東京23区で一番緑が多いのはどこ

特に東京に行ったこともないのですが、東京で一番緑が多いのはどこかがふと気になった次第です。将来引っ越しすることがあったら参考になるかもしれないので調べます。

やりかた

  1. 行政区域データを集める
  2. 土地利用データを集める
  3. 各行政区域について土地利用ごとの面積を算出
  4. 各行政区域について、自然っぽい土地利用の割合を算出

当然このようにすればできます。文字にすると簡単そうです。

データ収集

行政区域データは国土数値情報ダウンロードサイトから取得します。平成29年版の東京都のジオメトリデータで、ESRI Shapefile形式をダウンロードします。

土地利用データはALOSの高解像度土地利用土地被覆図の10mメッシュのgeotiffをダウンロードし、必要そうなところだけ展開します。

収集したデータは以下になります。geotiffだけで3GBくらいあります。

$ tree .
├── geotiff
│   ├── LC_N20E136.tif
│   ├── LC_N24E141.tif
│   ├── LC_N24E153.tif
│   ├── LC_N25E141.tif
│   ├── LC_N26E142.tif
│   ├── LC_N27E140.tif
│   ├── LC_N27E142.tif
│   ├── LC_N30E140.tif
│   ├── LC_N32E139.tif
│   ├── LC_N33E136.tif
│   ├── LC_N33E139.tif
│   ├── LC_N34E136.tif
│   ├── LC_N34E137.tif
│   ├── LC_N34E138.tif
│   ├── LC_N34E139.tif
│   ├── LC_N35E136.tif
│   ├── LC_N35E137.tif
│   ├── LC_N35E138.tif
│   ├── LC_N35E139.tif
│   ├── LC_N35E140.tif
│   └── ver1609_LC_GeoTiff.tgz
└── shapefiles
     ├── KS-META-N03-17_13_170101.xml
     ├── N03-17_13_170101.dbf
     ├── N03-17_13_170101.prj
     ├── N03-17_13_170101.shp
     ├── N03-17_13_170101.shx
     └── N03-17_13_170101.xml

DockerでPostGISサーバを立てる

行政区域の土地利用ごとの面積を計算するには、行政区域のジオメトリを土地利用ごとに切り出したジオメトリを用意する必要があります。 ジオメトリ演算といえばPostGISを使うようです。環境構築は面倒くさそうなので、dockerを使います。

以下のようなDockerfileを用意します。

FROM mdillon/postgis:latest
RUN rm -rf /docker-entrypoint-initdb.d/postgis.sh
EXPOSE 5432

/docker-entrypoint-initdb.d/にはコンテナ起動直後に実行するポスグレDB初期化スクリプトが配置されます。postgis.shはmdillon/postgisが用意したPostGIS初期設定スクリプトなのですが、不要なextensionとかもインストールするので削除します。 ちなみに/docker-entrypoint-initdb.d/内のスクリプトは名前順に実行されるっぽいので、適当な名前でオレオレスクリプトを配置すると実行順序が狂って死んだりするので注意が必要です。

$ docker build --no-cache -t postgis_calcgreen .
$ docker run \
  -e POSTGRES_DB=postgres \
  -e POSTGRES_USER=postgres \
  -e POSTGRES_PASSWORD=postgres \
  -v "$PWD":/usr/src/myapp \
  -w /usr/src/myapp \
  --name postgis_calcgreen \
  -p 5432:5432 \
  postgis_calcgreen

イメージつくってコンテナ起動します。パスワード等はテキトーです。/docker-entrypoint-initdb.d/のスクリプトがコンテナ起動直後に走るので、それを確認するために"-d"は付けないのが無難です。

動作確認します。psqlコマンド使います。postgresqlでは"\l"でデータベース一覧が取得できます。

$ export PGHOST=localhost
$ export PGDATABASE=postgres
$ export PGUSER=postgres
$ export PGPASSWORD=postgres
$ psql -c "\l"
                                 List of databases
   Name    |  Owner   | Encoding |  Collate   |   Ctype    |   Access privileges   
-----------+----------+----------+------------+------------+-----------------------
 postgres  | postgres | UTF8     | en_US.utf8 | en_US.utf8 | 
 template0 | postgres | UTF8     | en_US.utf8 | en_US.utf8 | =c/postgres          +
           |          |          |            |            | postgres=CTc/postgres
 template1 | postgres | UTF8     | en_US.utf8 | en_US.utf8 | =c/postgres          +
           |          |          |            |            | postgres=CTc/postgres
(3 rows)

ポスグレサーバが立ち上がりました。

データ投入とか

続いてpostgis extensionを追加します。

$ psql -c "create extension if not exists postgis;"
CREATE EXTENSION
$ psql -c "\dx;"
                                     List of installed extensions
  Name   | Version |   Schema   |                             Description                             
---------+---------+------------+---------------------------------------------------------------------
 plpgsql | 1.0     | pg_catalog | PL/pgSQL procedural language
 postgis | 2.3.2   | public     | PostGIS geometry, geography, and raster spatial types and functions
(2 rows)

できました。

行政区域シェープファイルはdockerコンテナ上でshp2pgsqlコマンドを使ってsqlに変換し、tmpという名のテーブルに登録します。-sはSRID、-iはInt4指定、-Iはgistインデックス作成、-Wは文字コード指定です。sqlに変換したら直接ローカルのpsqlコマンドに流して実行します。

tmpテーブルのカラムを調整し、全レコードを市区町村でgroup byしてST_Unionしてjapan_ksテーブルに登録します。tmpテーブルはお役御免とdropします。

$ docker exec -i $PGCONTAINER bash -c \
'shp2pgsql -s 4326 -D -i -I -W cp932 shapefiles/*.shp tmp' | psql >/dev/null
$ cat << EOS | psql
delete from tmp where n03_007 is NULL;

alter table tmp rename n03_001 to prefecture_name;
alter table tmp rename n03_007 to city_code;

alter table tmp add prefecture_code varchar(2);
update tmp set prefecture_code = substring(city_code,1,2);

alter table tmp add city_name varchar(60);
update tmp set city_name = coalesce(n03_003, '') || coalesce(n03_004, '');

create table japan_ks as
  select
    cast(prefecture_code as integer) as prefecture_code
    , prefecture_name
    , cast(city_code as integer) as city_code
    , city_name
    , sum(ST_Area(geography(geom))) as area
    , ST_Union(geom) as geom
  from tmp
  group by
    prefecture_code
    , prefecture_name
    , city_code
    , city_name
    ;
drop table tmp;
create index on japan_ks using gist (geom);
analyze japan_ks;
EOS

次にalos土地利用データを登録します。raster2pgsqlコマンドを使います。-Iはインデックス作成、-Cは「SRIDやピクセルサイズ等のラスタ制約を適用して、raster_columnsビューで適切な登録ができるようにします。 」らしいです。-Nはヌル値指定、-tはgeotiffをレコード単位で分割するときの、1レコード当たりのタイルサイズです。タイルサイズが大きいとジオメトリ演算で時間がかかりそうですが、小さいと登録に時間がかかります。

$ ALOS_GEOTIF=./geotiff/LC*.tif
$ docker exec -i $PGCONTAINER bash <<EOF | psql > /dev/null
raster2pgsql -I -C -s 4326 -N 0 -t 100x100 \
$ALOS_GEOTIF \
alos_rast
EOF

alos区分テーブルを下のsqlで作成します。

create table alos_code_name (
code integer primary key,
name varchar(40)
);
insert into alos_code_name (code, name) values
    (0,   '未分類'),
    (1,   '水域'),
    (2,   '都市'),
    (3,   '水田'),
    (4,   '畑地'),
    (5,   '草地'),
    (6,   '落葉広葉樹'),
    (7,   '落葉針葉樹'),
    (8,   '常緑広葉樹'),
    (9,   '常緑針葉樹'),
    (10,  '裸地'),
    (11,  '雪氷'),
    (253, 'その他'),
    (255, 'データなし')
;

alos_rastに格納したラスタデータをジオメトリに変換してalos_geoに格納します。行政区域ジオメトリと交差する土地利用ラスタのみ変換します。

create table alos_geo as
  select
    alos_geo.rid as rid
    , (alos_geo.geomval).val as code
    , (alos_geo.geomval).geom as geom
  from
    (select
      alos_rast.rid as rid
      , ST_DumpAsPolygons(alos_rast.rast) as geomval
    from
      (select ST_Union(geom) as geom from japan_ks) as whole_japan
      , alos_rast
    where
      ST_Intersects(whole_japan.geom, alos_rast.rast)
    ) as alos_geo
  ;
create index on alos_geo using gist (geom);
analyze alos_geo;

japan_ksとalos_geoのそれぞれのジオメトリのintersection(交差ジオメトリ)をjapan_alosテーブルに格納します。今回はalosのラスタをジオメトリに変換して(alos_geoテーブル)ST_Intersectionの引数に与えていますが、ST_Intersectionはラスタを直接引数にとることもできます(その場合は裏でラスタデータのST_DumpAsPolygonsが走る)。

create table japan_alos as
select
  japan_ks.city_code as city_code
  , japan_ks.city_name as city_name
  , alos_geo.code as alos_code
  , ST_Intersection(japan_ks.geom, alos_geo.geom) as geom
from
  japan_ks
  , alos_geo
where
  ST_Intersects(japan_ks.geom, alos_geo.geom)
;

結果

準備できました。可視化したり緑の多いところを調べたりします。Pythonを使います。GISの可視化といえばgeopandasなのだそうです。

%matplotlib inline
# ref. https://geohackweek.github.io/vector/06-geopandas-advanced/
import psycopg2

import matplotlib as mpl
import matplotlib.pyplot as plt
plt.style.use('bmh')

import pandas as pd
import geopandas as gpd

適当にmatplotlibの設定をします。

mpl.rc('font', family='Osaka') # Macだとこれで日本語が文字化けしないかもしれない
mpl.rc('legend', loc="lower left", fontsize='large')
# geopandas v0.3.0のlegend.locはbestでハードコードされているのでrcで指定しても無意味

PostGISサーバのアクセス情報です。

db_conn_dict = {
    "host":"localhost",
    "port":5432,
    "dbname":"postgres",
    "user":"postgres",
    "password":"postgres"
}

試しに行政区域ジオメトリを可視化します。

with psycopg2.connect(**db_conn_dict) as conn:
    japan_ks = gpd.read_postgis(
        "SELECT * FROM japan_ks", 
        conn,
        geom_col='geom',
        crs={'init': u'epsg:4326'}
    )

japan_ks.plot(column='city_code', cmap='tab20c', categorical=True, figsize=(14, 8))
plt.savefig("japan_ks.png", bbox_inches="tight", pad_inches=0.0, transparent=True

f:id:kazufusa1484:20171009220641p:plain

東京は広いですね。よくわかりません。島嶼を除きます。

japan_ks.plot(column='city_name', cmap='tab20c', categorical=True, figsize=(14, 8))
plt.ylim([35.45,35.95])
plt.xlim([138.85,140])

for idx, row in japan_ks.iterrows():
    plt.annotate(
        s=row.city_name,
        xy=row.geom.representative_point().coords[0],
        horizontalalignment='center'
    )

plt.savefig("japan_ks.zoom.png", bbox_inches="tight", pad_inches=0.0)

f:id:kazufusa1484:20171009220650p:plain

ラベルが少し見にくいですがそれらしいような気がします。なぜかbbox_inches="tight"が無効になっています。

土地利用別の可視化もしてみます。つまりjapan_alosテーブルを可視化します。下はカラーマップです。

alos_color_palette = {
    '00_未分類':     [0  /255 ,0  /255 ,0  /255 ],
    '01_水域':       [0  /255 ,0  /255 ,100/255 ],
    '02_都市':       [255/255 ,0  /255 ,0  /255 ],
    '03_水田':       [0  /255 ,128/255 ,255/255 ],
    '04_畑地':       [255/255 ,193/255 ,191/255 ],
    '05_草地':       [255/255 ,255/255 ,0  /255 ],
    '06_落葉広葉樹': [128/255 ,255/255 ,0  /255 ],
    '07_落葉針葉樹': [0  /255 ,255/255 ,128/255 ],
    '08_常緑広葉樹': [86 /255 ,172/255 ,0  /255 ],
    '09_常緑針葉樹': [130/255 ,172/255 ,117/255 ],
    '10_裸地':       [128/255 ,100/255 ,0  /255 ]
}

def alos_colors(alos_names):
    return mpl.colors.ListedColormap(
        [alos_color_palette[x] for x in alos_names],
        name="alos_colors",
        N=len(alos_names))
query = """
select
  japan_alos.geom as geom
  , japan_alos.city_name as city_name
  , alos_code_name.code as alos_code
  , to_char(alos_code_name.code, 'FM00') || '_' || alos_code_name.name as alos_name
from
  japan_alos
  left join alos_code_name on (japan_alos.alos_code = alos_code_name.code)
;
"""

with psycopg2.connect(**db_conn_dict) as conn:
    japan_alos = gpd.read_postgis(
        query, 
        conn,
        geom_col='geom',
        crs={'init': u'epsg:4326'}
    )
japan_alos.plot(
    column='alos_name',
    cmap=alos_colors(sorted(japan_alos.alos_name.unique())),
    categorical=True,
    legend=True,
    figsize=(14, 8)
)
plt.ylim([35.45,35.95])
plt.xlim([138.85,140.25])
plt.savefig("japan_alos.zoom.png", bbox_inches="tight", pad_inches=0.0)

f:id:kazufusa1484:20171009220635p:plain

ひとまず東京都できちんと切り出されていることは確認できました。

さて、市区町村ごとのalosの落葉広/落葉針/常緑広/常緑針葉樹の割合を算出します。

greenratio_query = """
select
    japan_alos.city_code
    , min(japan_alos.city_name) as city_name
    , sum(ST_Area(geography(japan_alos.geom))) / min(japan_ks.area) as green_ratio
from
    japan_alos
    left join alos_code_name on (japan_alos.alos_code = alos_code_name.code)
    left join japan_ks on (japan_alos.city_code = japan_ks.city_code)
where
    alos_code between 6 and 9
group by
    japan_alos.city_code
order by
    green_ratio desc, japan_alos.city_code
;"""

with psycopg2.connect(**db_conn_dict) as conn:
  greenratio = pd.io.sql.read_sql(greenratio_query, conn)

top10です。西多摩郡島嶼です。西多摩郡檜原村というところがほとんど森です。

greenratio.head(n=10)
city_code city_name green_ratio
0 13307 西多摩郡檜原村 0.992630
1 13308 西多摩郡奥多摩町 0.965666
2 13382 御蔵島村 0.915994
3 13361 大島町 0.836055
4 13362 利島村 0.832088
5 13305 西多摩郡日の出町 0.825971
6 13205 青梅市 0.763871
7 13364 神津島村 0.746610
8 13363 新島村 0.743949
9 13228 あきる野市 0.735391

可視化してみますが、やはりほとんど森です。

query = """
SELECT
    japan_alos.geom as geom
    , japan_alos.city_name as city_name
    , alos_code_name.code as alos_code
    , to_char(alos_code_name.code, 'FM00') || '_' || alos_code_name.name as alos_name
FROM
    japan_alos
    left join alos_code_name on (japan_alos.alos_code = alos_code_name.code)
where
    city_name = '西多摩郡檜原村'
;
"""

with psycopg2.connect(**db_conn_dict) as conn:
    hinoharamura = gpd.read_postgis(
        query, 
        conn,
        geom_col='geom',
        crs={'init': u'epsg:4326'}
    )

hinoharamura.plot(
    column='alos_name',
    cmap=alos_colors(sorted(hinoharamura.alos_name.unique())),
    categorical=True,
    legend=True,
    figsize=(14, 8)
)
plt.xlim([139,139.22])
# plt.gca().get_xaxis().get_major_formatter().set_useOffset(False)
plt.savefig("hinoharamura.png", bbox_inches="tight", pad_inches=0.0)

f:id:kazufusa1484:20171009220629p:plain

東京23区のランキングです。千代田区が1位ですね。

greenratio[greenratio.city_code < 13200].reset_index(drop=True)
city_code city_name green_ratio
0 13101 千代田区 0.110995
1 13113 渋谷区 0.092015
2 13103 港区 0.075061
3 13105 文京区 0.061869
4 13112 世田谷区 0.052596
5 13120 練馬区 0.046121
6 13115 杉並区 0.045913
7 13104 新宿区 0.037274
8 13109 品川区 0.035927
9 13110 目黒区 0.034893
10 13106 台東区 0.031348
11 13119 板橋区 0.028702
12 13108 江東区 0.024454
13 13117 北区 0.019062
14 13116 豊島区 0.016745
15 13111 大田区 0.016584
16 13102 中央区 0.016507
17 13114 中野区 0.016115
18 13122 葛飾区 0.015348
19 13123 江戸川区 0.013803
20 13121 足立区 0.008889
21 13107 墨田区 0.005948
22 13118 荒川区 0.003172

千代田区も同様に可視化してみました。緑を囲むお堀のような構造が見えます。これはつまり旧江戸城現皇居ですね。

f:id:kazufusa1484:20171009220555p:plain

皇居はずるいように思えたので、23区で2番目に緑が多い渋谷区です。こちらも緑色のオーバルが見られますが、これは代々木公園のようです。

f:id:kazufusa1484:20171009220709p:plain

というわけで23区で一番緑が多いのは千代田区、皇居は自由に入れないので除外すべきと思えば渋谷区となりました。

PostGISとかSQLは難しいので、上の手続きをgeopandasだけでできると大変素敵ですね。

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

階層ベイズモデルで勝敗データからプロ棋士の強さを推定する,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%も理解できない。