StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

「StanとRでベイズ統計モデリング」松浦健太郎 という本を書きました

Stan 書評

僕が筆者なので、この記事は書評ではなく紹介になります。まずこの本はRのシリーズの一冊にもかかわらずStanという統計モデリングのためのプログラミング言語の方がメインです。このようなわがままを許してくれた、ゆるいふところの深い石田先生と共立出版には感謝しかありません。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

目次と概要

共立出版のページを見てください。GitHubのリポジトリもあります。

前提とする知識

「はじめに」の部分で触れていますが、確率と統計の基本的な知識はある方、R（やPython）で簡単なデータ加工や作図が一通りできる方を想定しています。そのため、確率分布なんて聞いたことがない、プログラミングがはじめて、Rがはじめて、という方が読み進めるのは厳しいかもしれません。なお、Rの基本的な関数しか出てこないので、PyStanとmatplotlib（あるいはSeabornなど）でやるわっていうPythonユーザの方にも十分に読む価値があると思います。

Pythonユーザのための追記

この本を読んで習得できるもの

「統計モデリングの考え方」と「Stanの使い方」の二点です。

統計モデリングの考え方

基本的にデータ解析には「正解」がありません。検定の前提条件にせよ、統計モデルにせよ、機械学習のモデルにせよ、すべてのモデルは仮定にすぎないからです。しかし、ルール無用ではありません。得られた結果が有用であるため、統計モデリングにも沿うべき指針や考え方があります。統計モデリングは従来の検定ベースの統計と比べると、得意分野も考え方も大きく異なります。筆者は仕事がら普段は検定も使っていますが、「ある現象を理解したい・知識を獲得したい・予測したい」といったデータ解析の主目的については、検定では満足できない場合が多いです。その場合、統計モデリングがよい選択肢になります。この本では統計モデリングをきちんと習得できるように、考え方や手順といった正解がないものについても大胆に筆者の主張を書きました。

Stanの使い方

最近Stanを取り上げたデータ解析の本は増えてきています。僕が確認している範囲だけでも以下があります（刊行順）。

しかし、これらの多くは重回帰のような簡単な解析を数ページで紹介しているだけだったり、Stanの文法などは付録だったりしました。「Stanをしっかり学びたい。より美しく、より速く動くように書きたい。」という人にとって十分に満足できる本はこれまでなかったと言えます。しかし、本書はこのニーズにほぼ100%応えた、日本語で読めるはじめての本格的なStan解説書です。

このブログとの関係

いくつかの記事については本の後半の例題とオーバーラップがありますが、基本的には書き下ろしです。本では用語の説明や結果の解釈などの基礎から順番に説明しています。このブログは僕のメモを兼ねているのでやや難しい題材やStanコードが多いですが、本ではそんなことありません。

その他の特徴

例題からはじまる

なるべく具体的なイメージをもって分かりやすくなるように、多くの解説は例題データの解析からはじまります。

図が多い

統計モデリングを含むデータ解析は可視化と密接に関係しているので、図をふんだんに載せています。図と作図のRコードもGitHubにほぼすべてありますので、よろしければ見てください。ただし、作図は各自が好きなパッケージを使えばよいと思っているので作図のRコードの解説は割愛しました。すみません。

数式は難しくない

ベイズ統計の本では1ページまるまる難解な式変形という本も少なくないです。しかし、この本はソースコードとその解説が多めですが、数式は少なめです。式変形も理系の大学一二年生なら引っかかることなく追うことができるでしょう。

Stanの最新版に対応していく

書籍は2.11対応ですが現時点の2.15にもほぼ対応しています。今後もバージョンアップしたら、GitHub上でどの記述が古くなったか書きます。またソースコード自体もアップデートする予定です。

他人からの書評や感想 (追記)

献本勢が多いので褒める言葉は話半分でよいと思いますが、他人から見てよかった章の紹介などは参考になると思います。

この他にも色々コメントいただいております。ありがとうございます。

読書会・勉強会

読書会が開かれることになりました。大変ありがたいです。

なお、愛称は「アヒル本」になった模様です（表紙のブロックはアヒルなのです）。

各章の紹介

以下では実物を手に取って見ることができない人のために、どんなことが書かれているか簡単に紹介したいと思います。

はじめに

本書のあらすじが載っています。

f:id:StatModeling:20201106183153p:plain

1～3章は理論編です。数式が多いわけではなく用語の確認の意味合いが強いです。4～5章はStanの入門編です。6章以降は発展編です。6・8・11・12章が本流で、モデルのレパートリーを増やす章になっています。対して、7・9・10章はモデルを改善する章となっています。特に発展編では読者はStanの強力さを実感できると思います。

1章 統計モデリングとStanの概要

「統計モデリングとは何なのか？目的は？メリットは？」という点について筆者なりの考えを述べました。またStanなどの確率的プログラミング言語を使うメリットを簡潔に説明しています。一部抜粋します。

確率的プログラミング言語とは「様々な確率分布の関数や尤度の計算に特化した関数が豊富に用意されており、確率モデルをデータにあてはめることを主な目的としたプログラミング言語」である。ユーザーはモデルをプログラミングコードで記述し、データを渡すだけでよい。すると確率的プログラミング言語の方でほぼ自動的にパラメータの値を推定してくれる。このようにモデルの記述と難しい推定計算を分離することによって、モデルの可読性が上がり、バグの混入が激減し、解析者はモデルの試行錯誤に専念できるようになる。特に多数のモデルを試行錯誤する状況で確率的プログラミング言語は真価を発揮するのである。

2章 ベイズ推定の復習

確率分布や「 y \sim Poisson(\lambda)」などの基本的な用語や記法の確認からはじまります。そのあとでベイズ統計とMCMCに関する用語を簡潔に説明しました。扱った用語は以下の通りです。

3章 統計モデリングをはじめる前に

まず一般にデータ解析に必要となる前準備について説明しています。そのあとで統計モデリングの手順を筆者なりに以下のように定型化しました。

  • 解析の目的
  • データの分布の確認
  • メカニズムの想像
  • モデル式の記述
  • Rでシミュレーション
  • Stanで実装
  • 推定結果の解釈
  • 図によるモデルのチェック

本書の以降の例題はなるべくこの手順に沿って解析をすすめています。もちろん、現実のデータ解析はこの手順に尽きているわけではなく、あくまでも最低限踏むべきステップの目安です。また、事前知識の役割や、誤解の多いモデルの「正しさ」についても少し触れました。

4章 StanとRStanをはじめよう

インストール方法、基本的な文法、targetとlp__を説明した後で、じっくり単回帰の問題を扱います。この章以降ではStanコード・Rコードとそれらの説明という部分が多くなります。推定結果の見方、収束の判断、trace plotの見方、{ggmcmc}パッケージの使い方、MCMC設定の変更、MCMCサンプルの使い方を説明しています。

5章 基本的な回帰とモデルのチェック

重回帰・二項ロジスティック回帰・ロジスティック回帰・ポアソン回帰を扱います。また図によるモデルのチェック方法を数通り紹介しています。使用した図をいくつか載せておきます。詳しくは本を読んでください。

f:id:StatModeling:20201106183157p:plain

f:id:StatModeling:20201106183201p:plain

f:id:StatModeling:20201106183206p:plain

6章 統計モデリングの視点から確率分布の紹介

確率分布の軽い紹介はどんな本にも載っていて、常々冗長だ退屈だと感じていました。しかし、統計モデリングにおいて個々の確率分布は必要不可欠なパーツです。そこで、従来の本にはないような「統計モデリングの視点から」の部分になるべく注力して16個の確率分布を紹介しました。

7章 回帰分析の悩みどころ

交互作用、対数をとるか否か、非線形の関係、打ち切り、外れ値など、ふつうの統計の教科書にはあまり載っていませんが、実際のデータ解析では必ずといっていいほど直面する悩みどころを取り上げました。

8章 階層モデル

階層モデルはグループ差や個人差をうまく扱うための一手法です。ゆっくり導入をした後で応用方法になじめるように、複数の階層を持つ場合や、非線形モデルやロジスティック回帰の階層モデルを取り上げました。

9章 一歩進んだ文法

Stanをより美しく、より速く動くように必要な一歩進んだ文法を解説しました。Stanで用意されている型（かた）の説明やベクトル化、行列演算、パラメータの制約、欠測値がある場合などを扱っています。

10章 収束しない場合の対処法

統計モデリングの最大の難所は、MCMCが収束しないことだと言っても過言ではありません。しかしながら、その対策を系統的に扱った書籍や情報はあまりありません。ここではその対策を大きく4つの節に分けて説明しています。個人的な経験からは、10.1節の識別性の問題（多重共線性を含む）と10.2節の弱情報事前分布で収束しない場合の8割は解決できそうです。

11章 離散値をとるパラメータを使う

Stanの最大の弱点は離散値をとるパラメータが直接的に使えないことです。この章では、log_sum_exp関数と周辺化消去（marginalizing out, summing out）で解決する方法を説明します。応用例として、混合正規分布、ゼロ過剰ポアソン分布、Latent Dirichlet Allocationを扱います。変分ベイズ法の一実装であるADVIも使ってみます。

12章 時間や空間を扱うモデル

時系列データを解析するにあたり、応用範囲が広く、解釈がしやすく、拡張性が高い「状態空間モデル」を取り上げました。特に日次の売り上げデータのように、時刻が離散的で等間隔なデータの場合に使いやすいモデルです。季節調整項や変化点検出なども扱います。後半では、時間構造と空間構造の等価性を説明し、マルコフ場モデル（CARモデルとほぼ同じもの）を使った例題を扱っています。

データ解析で割安mobile PCを探す

現実のデータ Stan

この記事の続編です。一緒にやろうという人がなかなか現れないので、一人でたたき台を作りました。

目的

目的は機能の割にお得な割安mobile PCを探すことです。mobile PCの厳密な定義はないのですが、ここではディスプレイが12型～14型で重さが1kg前後としました。また、各社の最新モデルだけを対象としました。

データの取得方法

メーカーを決める→本気で買うつもりで公式サイトと価格comを比較して安い方にする→人力スクレイピング です。現時点では公式サイトも多く、スクレイピングのコードを書いても労力のもとが取れないので10時間ほどかけて人力スクレイピングして集めました。

データの内容

8社・10モデルで44商品です。おすすめモデルを中心にしました。生データを置いておきます（これとこれ）。

次の統計モデリングで使用する、PCの機能を表す説明変数は18個考慮しました。CPU・メモリ・SSD・PCIe・ディスプレイのサイズ/解像度/光沢非光沢・重さ・バッテリ持続時間・キーボードバックライト・LANポート・SIMポート・USB3.0の数・USB3.1の数・各種出力などです。前処理はそこそこ必要で、例えばCPUはこのサイトのベンチマークスコアの値に変換して使いました。また、キーボードの打ちやすさ、画面の見やすさ、公称でないバッテリー持続時間などはこちらのサイトから取得して、一部を説明変数として使いました。目的変数は税込み価格です。

簡単な可視化は記事の最後にある発表資料を見てください。

統計モデル

切片・説明変数の項・ブランド（会社名の影響）の項からなるシンプルな回帰です。ブランドは正規分布に従うとし、その標準偏差は弱情報事前分布を設定し、データから推定しました（階層モデル）。例えば「Think Padのキーボードが打ちやすい」という特徴があっても、今回は会社の影響と切り分けできないのでブランドの影響に組み込まれることになります。

計算方法

StanとRで計算しました。記事の最後の方にコードを載せておきます。計算時間はSurface Pro 3で1chainあたり約20秒ほどです。

結果

説明変数の影響

f:id:StatModeling:20201106182810p:plain

黒い点が推定した事後分布（からのMCMCサンプル）の中央値で、横に伸びている線が95%ベイズ信頼区間です。

影響があると断言できそうなのは、CPUとメモリ容量とSSD容量ぐらいでした。重さは軽い方が高い傾向はあるのですが、95%区間がゼロをまたいでおり、そうも断言できない結果となりました。

ブランドの影響

f:id:StatModeling:20201106182814p:plain

凡例は「説明変数の影響」と同じです。大差なしという結果です。あると思っていたのでやや残念です。

実際の価格 vs. 「潜在的な価値」

f:id:StatModeling:20201106182819p:plain

凡例は「説明変数の影響」と同じです。直線y=xの上側が割安な商品、直線y=xの下側が割高な商品です。ほとんどの商品が直線y=xの近くにあるので、目立つ割安・割高な商品はないと言えます。市場原理はなかなか強力のようです。どれも値段相応の価値があると言えるので、好きなもん買えばいいと思います。

しかしよく見ると、富士通の商品はMCMCサンプルの中央値がすべて割安側にあり、また（性能に関係ない）ブランド効果も比較的低いので、会社として見ると割安な商品を出す傾向にあると思います。たしかに公式サイトの訳あり商品はなんかいつも安い気がする。

割引%offのBest3・Worst3

目立った割安・割高な商品はないですが、例えばMCMCサンプルの中央値ベースで割引%offのBest3・Worst3なんかも求めることができます。

best1best2best3worst1worst2worst3
price103058156800109938219427139655167184
潜在的な価値12.2万円17.40万円12.09万円18.76万円12.33万円15.10万円
%off15.59.99.1-16.9-13.2-10.7
search_date201607172016071720160717201607172016071720160717
search_sitekakakukakakukakakukakakuofficialofficial
urlhttp://kakaku.com/item/K0000872152/http://kakaku.com/item/K0000855827/http://kakaku.com/item/K0000752290/http://kakaku.com/item/K0000855828/http://shopap.lenovo.com/jp/notebooks/thinkpad/x-series/x1-carbon/#tab-customizehttp://www.apple.com/jp/shop/buy-mac/macbook-air
company_nameApplePanasonicApplePanasonicLenovoApple
model_nameMacBookAirLet'snoteSZ5MacBookAirLet'snoteSZ5ThinkPadX1CarbonMacBookAir
CPUCorei5-5250UCorei5-6300UCorei5-5250UCorei5-6300UCorei5-6200UCorei7-5650U
memory844848
SSD128128256256128256
PCIe101001
display_size13.312.113.312.11413.3
display_res_width144019201440192019201440
display_res_height900120090012001080900
display_touch000000
display_glare101001
display_vote000010
keyboard_light101011
keyboard_vote000010
battery_size54NA54NA5254
battery_time11211.512139.812
battery_time2NANANANA5.62NA
weight1.350.8751.350.8491.181.35
DVD000000
SIM010000
LAN_port010100
USB30333333
USB31000000
output_HDMI111111
output_VGA010100
SDcard111111
MSoffice000000

best1とbest3はMac Book Airのローエンドな商品です。これらは少々古い割に人気があるので、取り扱っている店が多く、価格comで激しい競争にさらされてお買得になっていると思われます。一方で、Appleの公式サイトで選んだハイエンドな商品はworst3になっています。Appleの公式サイトで買うのは割高の傾向があると言えるでしょう。

best2とworst1はLet's note SZ5です。ともに価格comで扱っている店も十分にあるのにこの差があるところが面白いです。1商品にしぼって価格comを信じるのではなく、データ解析の関係者ならメタな視点でお買い得な商品を選びたいものです（まあ統計モデルもお買い得かも仮説にすぎないのですが）。

Future work

価格comのWeb APIやAmazon APIを使って、商品を最新モデルに限らず、さらに扱っている店舗の影響や価格推移などのデータを集めるのは考えられます。ただし、取り扱っている店の少なさ、限定モデルかどうか、などをチェックしてフィルターする必要がありそうです。しかし、データが増えれば、「ある会社の製品は早く値段が落ちる」なんてことも分かりそうな気がします。データを定期的に取得するのが面倒なんだよなぁ。

ソースコード

Stanコードは以下です。

data {
  int N;
  int D;
  int K;
  matrix[N,D] X;
  vector[N] Y;
  int N2K[N];
}

parameters {
  vector[D] beta;
  vector[K] brand;
  real<lower=0> s_Y;
  real<lower=0> s_brand;
}

transformed parameters {
  vector[N] mu;
  vector[N] mu_plus_brand;
  mu = X*beta;
  mu_plus_brand = mu + brand[N2K];
}

model {
  beta[1] ~ student_t(4, 0, 3);
  beta[2:(D-1)] ~ student_t(4, 0, 1);
  s_brand ~ student_t(4, 0, 1);
  brand ~ normal(0, s_brand);
  Y ~ normal(mu_plus_brand, s_Y);
}

実行するRコードは以下です。

library(rstan)

d <- read.csv(file='input/data.csv', stringsAsFactors=FALSE)
d_cpu <- read.csv(file='input/cpu_info.csv', skip=1, stringsAsFactors=FALSE)
rownames(d_cpu) <- d_cpu$CPU
company_names <- c('Lenovo', 'Dell', 'NEC', 'Panasonic', 'Fujitsu', 'VAIO', 'HP', 'Apple')
d$company_name <- factor(d$company_name, levels=company_names)
d$cpubenchmark <- d_cpu[d$CPU, ]$cpubenchmark/4000
d$price <- d$price/100000
d$display_res <- (d$display_res_width * d$display_res_height)/1920/1080
d$battery_time1 <- d$battery_time1/20
d$memory <- d$memory/16
d$SSD <- d$SSD/512
d$ports <- d$LAN_port

usb_cvt <- c('0'=0, '1'=0.5, '2'=0.8, '3'=1)
d$USB30_val <- usb_cvt[as.character(d$USB30)]
d$USB31_val <- usb_cvt[as.character(d$USB31)]

Y <- d$price
gr <- d$company_name
use_cols <- c('cpubenchmark', 'memory', 'SSD', 'PCIe', 'display_size', 'display_res', 'display_touch', 'display_glare', 'keyboard_light', 'battery_time1', 'weight', 'DVD', 'SIM', 'ports', 'USB30_val', 'USB31_val', 'output_HDMI', 'MSoffice')
X <- data.frame(1, d[ , use_cols])
data <- list(N=nrow(X), D=ncol(X), K=nlevels(gr), X=X, Y=Y, N2K=as.numeric(gr))

stanmodel <- stan_model(file='model/model.stan')
fit <- sampling(stanmodel, data=data, seed=1)

発表資料

9/7にYahoo!でこの内容でLTしました。資料は以下です。

「Python機械学習プログラミング」 Sebastian Raschka（著）, 株式会社クイープ（訳）, 福島真太朗（監訳）

書評 Python

僕はベイズ統計モデリングをはじめる前（5年ほど前）までは主に機械学習をしていました。その頃は平易な成書はあまりなくて、サポートベクターマシンの理論の難しい本を読んだり、Weka本（当時はこれ）を読みながら実装していたことを思い出します。PythonでもSVM-RFEを書いたりしてました。しかし、時は流れ、Pythonからscikit-learnという機械学習用ライブラリや深層学習を手軽に使うことができるようになり、気づいたらPythonは機械学習に必要不可欠な言語になっていました。この本はそんな機械学習に特化したPythonの使い方を理論と実装の両面から平易に丁寧に説明しています。理論は理系学部生なら理解できるぐらいで、実装はPythonやnumpyを少し触ったことがある人なら分かるぐらいのレベルです。いつの間にかこのような読みやすい機械学習の和書が出ているのは感慨深いです。

Python機械学習プログラミング 達人データサイエンティストによる理論と実践 (impress top gear)

Python機械学習プログラミング 達人データサイエンティストによる理論と実践 (impress top gear)

この本のおおまかな内容は、教師あり学習/教師なし学習の概論（1章）、分類（ロジスティック回帰/サポートベクトルマシン/決定木）（3章）、前処理（4章）、次元削減（PCAなど）（5章）、モデル評価とハイパーパラメータのチューニング（6章）、アンサンブル学習（7章）、応用例（8章・9章）、回帰分析（10章）、クラスタ分析（11章）、ニューラルネットワーク（12章・13章）です。原著に対応したコードがGitHubにあります。書籍中のコードはコメント部分も日本語になっているので読みやすいです。

以降ではよかった点について簡単に触れたいと思います。

理論と実装のバランス

実装の方に重きがありますが、多くのトピックで実装に移る前に簡潔な理論の説明があります。バランスが絶妙です。原著の著者はミシガン州立大学の博士課程の方のようで、時に一見ギャップのある式変形やサイエンスでは常識となっている単語（要素還元主義とか）がいきなり出てきますが、監注（監訳の注）でそういうギャップはこれでもかってほど埋められています。

scikit-learnに詳しくなれる

これまでにscikit-learnについて詳しくなろうと思ったら、公式のドキュメントを英語で読むか、日本語なら無限猿さんのブログ記事しかまとまった情報がありませんでした。しかし、この本はざっと200ページぐらいにわたってscikit-learnを使っています。

派手な分類やクラスタリングと異なり、4章・6章・7章のような一見地味だけど重要な情報は類書や日本語で書かれたWebページには詳しい説明がない場合が多く、価値が高いと思います。

ゼロベースの実装が教育的

3章以降では効率的で汎用的な機械学習ライブラリであるscikit-learnを基本的に使います。しかし、現実の問題に取り組んでいるユーザはscikit-learnだけで閉じれるわけではありません。論文が出たばかりのアルゴリズムを試す必要や、scikit-learnでは手の届かないようなフローの追加をする必要に日々迫られます。しかし、この本ではPythonを使ってゼロベースでアルゴリズムを実装することもしばしば取り上げられており、非常に教育的です。

例えば、2章ではゼロから確率的勾配降下法を実装しています。5章ではゼロからカーネル主成分分析を実装しています。12章ではゼロからニューラルネットワークを実装しています。ちょっと凝ったアルゴリズムがこんなにシンプルに書くことができるのかと驚きます。

matplotlibに詳しくなれる

データサイエンスは可視化と切っても切れない関係です。機械学習も例外ではないと思います。この書籍ではとにかく結果をmatplotlibで可視化しており、グラフが多くて見てて楽しいです。若干可視化のコードの分量が多すぎる気もしますが、matplotlibのプロの使い方を学べるのでよしとしましょう。

補足

監訳の福島さんがブック・インサイド―『Python機械学習プログラミング』学び方ガイドを書いています。「本書の前提知識」「読み方のプラン」「参考文献」などが載っています。「読み方のプラン」とか親切だなぁ。本は分からんところは飛ばしつつも全部読むもんだと思ってましたわ。