2021-02-26

【書評】意思決定分析と予測の活用　基礎理論からPython実装まで

書評

本記事について

著者の「Logics of Blue」馬場さんより、2021/02/25発売の「意思決定分析と予測の活用　基礎理論からPython実装まで」をご恵送いただきました。ありがとうございます！

www.kspub.co.jp

経営工学やオペレーションズ・リサーチが専門の同僚がいるのですが、畑違いのためどういう専門領域なのかわかっておらず、理解したいなと思っていたところでした。本書の内容はまさにこれらの領域や、システム工学、決定科学などの領域で活用されるようです（本書p4より）。

初学者として拝読しましたが、非常に面白かったです。本書で取り上げられている例の一つに、「ベネフィットとコストを勘案すると、機械を何台稼働させればいいか」という意思決定問題があります。僕だったら漠然と「うーん、とりあえず、3台くらい動かしてみてから考えよっか」と決めたくなってしまいます。

でも好況不況や、それぞれの機械のランニングコスト、社会的需要などについての情報があれば、ある基準を最適化するような予測が導ける、ということが本書で解説されています。

考えてみればそりゃあそうで、僕自身の日常生活では意思決定のミスが起こっても、そこで生まれる金銭的損失はそんなに大きくありません。うっかり終電を逃してしまって、仕方ないからタクシー使うか、てへへ、くらいのもんです。

でも金融とか経営とか、動くお金が膨大になるような場面では、ほんの小さな選択肢の違いが、とんでもない利得／損失の違いに帰結することもありえますよね。そういった場面では感覚的に意思決定するわけにはいかないので、手元にある情報から計算して、ある基準を最適化する予測を導く必要があるということだと理解しました。

【追記】

「予測を導く」っていうのは、本書の内容を表す正確な表現じゃなかったかもしれません。色んな予測結果が得られているうえで、予測結果をどう活用するか、ということが本書の焦点のようです。

以下、簡単な書評です。

本書の特長

1.　簡単な事例で、豊富な理論を紹介

本書で解説されている内容は、基本的に「利得行列」と呼ばれる情報の一覧を手掛かりにしています。現実的な意思決定では、この利得行列が複雑になったり、状況を想像することが難しい場合もあるのだと思いますが、本書では非常にシンプルな事例を用いているため、本質的な部分が分かりやすいと思います。

本当に多くの理論や数式が紹介されるので、「この添え字jって何だっけ...」など混乱することもあると思いますが、事例がシンプルなのですぐに立ち戻って思い出すことができます。

対して、数式や理論は多く紹介され、かなり濃密な内容だという印象を受けました。例えば第2部のエントロピーの話などは難しい内容だと思います。この辺の内容は、（馬場さんの本ではありませんが）以下の書籍なども参考にするといいかもしれません。

www.asakura.co.jp

第3部になると、さらに数式の難易度が上がります。数式だけを見ると、「ｳｯ...」と気後れしてしまう人もいるかもしれませんが、ちゃんと地の文を読むと、丁寧な解説が行われています。「この数式のこの部分、どういう意味なんだ...」と思ったらその近くには「これは～～という意味です」と解説がついています。数式パートは紙と鉛筆を用意して理解し、後述するプログラミングパートではPCで実装する、という学び分けが良いのではないでしょうか。

馬場さんは、他の書籍でもそうですし、実際に講演を拝聴したこともあるのですが、本当に読者／聴衆を置いてきぼりにすることがありません。躓きそうな石があったらちゃんと除いてくれる、そんな気遣いを感じます。

2.　Pythonコード

馬場さんの本といえば、事例、数式、そしてプログラミングコードの3点セットでお馴染みです。恐らく本書で紹介されている全ての計算において、Pythonコードが掲載されていると思います。実際に手を動かしながら学べるので、自分が今何をやっているのか、理解しやすいと思います。

多くのライブラリを使いこなして計算するのではなく、numpyとpandasを組み合わせて、時には関数を自作して、数式を実装していく感じです。馬場さんの著書「Pythonで学ぶあたらしい統計学の教科書」の進め方に近いですね。

www.shoeisha.co.jp

注意が必要なのは、本書はPythonそのものの解説書ではないということです。初めてPythonを使う人を想定して、関数の定義方法や、numpyのndarray、pandasのDataFrameなどの簡単な説明はありますが、必要最低限にとどめられています。

また、本書はコードをJupyter Notebook上で実行することを想定していますが、Jupyter Notebook自体の説明はありません。

ただし、これらは説明が不十分と捉えるべきではないと思います。本書があくまで伝えたいのは理論の部分であり、それらを丁寧に解説した結果すでに400ページ近くの分量になっています。ここでPython自体の説明も紙面を割いてしまうと、本当に伝えたいところを削らざるを得ないという判断だと思います。

サポートページはしっかりと用意され、そのURLは本書内に記載されていますし、有益な参考図書も多く紹介されているので、全く心配はいりません。

logics-of-blue.com

Python自体を初めて触る人は、適宜、他の参考図書なども併用しながら読むと良いと思います。

まとめ

僕自身、この書籍がターゲットとする内容が専門ではないので、踏み込んだ書評ができませんでしたが、初学者目線では非常に勉強になりました。第4部では、選好や効用関数など、経済学でお馴染みの内容が紹介されるので、経済学に関心がある人は、この辺りの内容も役に立つかもしれません。

僕自身、これから何度もこの書籍を読み返し、勉強していきたいと思います。

Enjoy！

2021-01-26

【書評】データ分析のためのデータ可視化入門

書評

本記事について

訳者の方々から、2021/01/26発売の、「実践Data Scienceシリーズ　データ分析のためのデータ可視化入門」をご恵送いただきました。ありがとうございます！

www.kspub.co.jp

本書は、「Data Visualization A practical introduction」の邦訳です。

socviz.co

原著は過去に私費で購入していたのですが（写真左）、自分で読むのも人に貸すのも、和書のほうが圧倒的に楽なので、たいへん嬉しいです。

f:id:das_Kino:20210126094606p:plain

以下、書評です。

本書の特徴

1.　「可視化という作業」自体の解説

タイトルからも分かるように、本書はRのggplot2パッケージを用いた可視化の入門書です。この手のソースは、事例集やハンズオン形式のコード紹介に注力したものが多いですが、本書ではまず最初に、なんと40ページ近くも費やして、「可視化って奥が深いんだぞ！」ということを熱弁しています。

可視化は統計的解析とセットで用いられることが多いと思いますが、可視化を怠ると、統計的解析において深刻な問題が発生しえます。有名なのは、「記述統計量（例：平均、分散）が同じだが、データの分布がまるで異なることがある」という事例ですね。

journals.plos.org

本書でも「アンスコムの四重奏」という事例から始まり、様々な観点で可視化の重要性を訴えています。

個人的に面白いと思ったのは、「ポップアウト」や「ゲシュタルト」など、知覚・認知心理学的な観点でも説明している点です。情報の取捨選択に係る能力を注意（attention）と呼びますが、「頑張って注意を払わないと情報が取得できないグラフ」もあれば、「注意をさほど払わなくても情報を選択できるグラフ」もあります。ヒトのことを理解したうえで、ヒトが効率的に情報を獲得できるように工夫しよう、というメッセージが込められていて、第1章はとても読みごたえがあります。

2.　丁寧な事例解説

ggplot2は、ひとたびクセを掴むと、かなりスラスラ書けます。それは可視化の文法（Grammar of Graphics）に基づいて設計されているからです。しかし、細かい調整をしようとすると、めったに使わない関数や記法が必要となり、とたんに難しくなります。

本書では、基本的なggplot2の使い方を丁寧に解説しているだけでなく、けっこうマイナーな書き方まで紹介しています。マイナーというのは、「そんなんいつ使うねん」という意味ではなく、「ちゃんと説明しようとすると難しいから、多くの書籍では恐らく回避されていたところ」という意味です。

原著の「A practical introduction」という副題の通り、実践的な用途において必要となる、痒いところに手の届く書き方まで詳しく解説されています。例えば第6章では「予測の図示」という方法が紹介されています。散布図に対して近似直線を引くだけなら簡単なのですが、予測というのは直線的関係だけに限りませんし、実用的には推定結果も点推定だけでなく区間推定まで行いたいものです。そういったところまでちゃんと解説されているのが嬉しいです。

3.　地理データの描画　

これは訳者の一人、瓜生さんが何より得意とされている内容だと思います。ggplot2では、地理データの描画が可能です。例えば日本地図を、人口密度に応じて都道府県ごとに色を塗り分ける、などが可能です。最新のバージョンでは、geom_sf()という関数が存在し、地理空間システム（GIS）データを扱いやすくなっています。

uribo.hatenablog.com

恐らく原著が出たタイミングでは、まだgeom_sf()が実装されていなかったのだと思うので、本書では既存のgeom_polygon()を用いて地理データが描画されています。しかし地理データの色を塗り分けたり、地理的な空間配置を考慮したグラフを描いたりする際に、本書の内容は最新のggplot2でも活用できると思います。

個人的には、地理データの描画はすごく苦手だったので、こうして体系的な解説が日本語で読めるというのは、非常にありがたいです。

4.　訳者の技量

僕は（もう閉会しましたが）Hijiyama.Rなど、西日本で勉強会に参加することが多く、Tokyo.Rなど東日本の勉強会に参加したことがありません。そのため実は、本書の訳者のお三方とは、直接の面識はないのですが、皆さん様々な解説記事・資料をインターネット上に公開されているので、日頃から「なんてすごい人たちなんだ」と思っていました。

そんな方々が翻訳されているので、日本語としての自然さや、訳語の適切さについては、言うまでもなく優れていると思いました。ただ何より僕が感動したのは、「翻訳者の判断で、原著とは違うコードを使用したこと」です。

tidyverseパッケージ群は、開発が盛んであるがゆえに、関数がコロコロ変わりがちです。335ページの脚注を見ると、「原著では（dplyrパッケージの）sample_n()関数を使用していたが、dplyrバージョン1.0.0から追加されたslice_sample()関数に書き換えた」とあります。最新の仕様をフォローしているからこそ出来ることですよね。

まとめ

上記のように、本書は単なる「ggplot2の解説書」ではなく、「ggplot2を用いた、可視化の解説書」になっています。識者の方々が翻訳を担当されたことにより、最新の仕様も考慮した訳書になっています。

本書はかなり詳しく書かれているぶん、上記のように深いところまで説明しているので、もしかしたら「うわああ、覚えることが多い～～」と焦ってしまう人もいるかもしれません。まず本書の第1章を読んで可視化の重要性について理解を深めたのち、いったん別の情報源でggplot2を簡単に学び、そのあとで本書に戻ってくるのでもよいかなと思いました。

Hadley Wickham著の「R for Data Science」はお勧めです。

r4ds.had.co.nz

あと、書評に拙著を紛れ込ませるのは下品だと思うのですが（本当にごめんなさい）、「RユーザのためのRStudio［実践］入門」も入門書としてお勧めです。現在、内容の刷新＋内容の大幅追加を行った第2版の改訂作業を行っている最中ですので、遠からずご案内できると思います！

gihyo.jp

すいません、最後に宣伝をしてしまいましたが、「実践Data Scienceシリーズ　データ分析のためのデータ可視化入門」は本当にお勧めです。皆さんも是非！　Enjoy！

2021-01-15

【書評】「Rで学ぶ統計的データ解析」

書評

本記事について

講談社サイエンティフィク様より、データサイエンス入門シリーズ「Rで学ぶ統計的データ解析」を教科書見本としてご恵送いただきました。

www.kspub.co.jp

実はデータサイエンス入門シリーズはほぼ全巻購入していて、大好きなシリーズです。

www.kspub.co.jp

書評というと偉そうですが、この本も本当に良い本だなと思ったので、簡単に感想を書きたいと思います。

（フルカラーで3000円って破格ですよね...）

本書の特徴

1.　最小限のR環境

最近はもう、RはRStudio上で実行するのが標準的だと思いますし、tidyverseパッケージ群（例えばデータフレーム操作はdplyr、可視化はggplot2）を活用することが増えてきたと思います。

個人的にはRStudio＋tidyverseパッケージ群が好きですが、導入はそれなりに面倒くさいとも思います。RStudioは多機能が故にどこを見ればいいか最初は分かりづらいし、tidyverseパッケージ群はインストールに失敗することもあります。あと開発が盛んなので、関数が増えたり減ったり変わったりもしますし。

一方本書では、RStudioは使わず素R環境で説明が行われていますし、特別なパッケージのインストールもほとんど求めていません。実はRは、デフォルトでもかなり多くの計算用関数・可視化用関数があるんですよね。第1章・第2章で、baseのRだけで出来ることを、かなり詳しく紹介してくださっています。apply()ファミリーなど、ちょっと癖があるけど使いこなしたら便利な関数もちゃんと説明されているのが嬉しいです。

特に第2章の「データの可視化と要約」のところは、かなり勉強になりました。僕はRを本格的に使い始めた直後からggplot2を勉強したので、実はbaseのRで可視化する方法をほとんど理解していなかったんですが、「こんなに色んなことできたんだ...」と驚きました。

もちろんパッケージを全く利用していないという意味ではなくて、例えばクラスター分析でclustrdパッケージを使うなど、要所要所で必要なパッケージの紹介もなされています。

2.　明確な道筋

データ分析の教科書には、色々な分析方法が紹介されていますが、初学者のうちはそれぞれの分析がどのような関係にあるのかを理解するのは容易ではないと思います。例えば「回帰分析」と「分散分析」は多くの教科書に載っていると思いますが、実際には互いに密接な関係があるにもかかわらず、これらの関係はあまり明確に言及されていないことも多いと思います。

一方本書では、

第3章：回帰分析
第4章：回帰分析（特に、正則化法を用いたリッジ回帰など発展的方法）
第5章：判別分析（被説明変数が質的変数の場合における回帰分析のようなもの）
第6章：ロジスティック回帰モデル（判別分析と関係した分析）
第7章：決定木に基づく判別モデル（判別分析やロジスティック回帰モデルとは違う強みがあることを紹介）
第8章：主成分分析（判別分析と似ているところはあるが、情報の圧縮という、別の目的のため）
第9章：クラスター分析（判別分析などと同様、分類が可能だが、教師あり／なしという違いがある）

という流れで、明確に筋の通った話が一貫して続いています。各章で全然違う分析を紹介しているという印象は全くなく、分析同士がどのような関係にあるのかが非常に明確に説明されていると思います。

また第10章以降は、ブートストラップ法やシミュレーションの実演が行われており、「こういう魔法みたいな手法／法則があるんだよ」ではなく、その挙動を実際に確かめられるようになっています。

3.　「まずは実行しよう、数理はそれからだ」という帯

帯に書いてあるコメントです。これだけ読むと、Rのマニュアルに特化した本なのかと思ってしまうかもしれませんが、全くそんなことはありません。プログラミング一辺倒でもなく、数式一辺倒でもなく、バランスよく両者が併記されていると思います。

バランスよくと書きましたが、これは分量の話です。

僕自身は数式をざっと眺めて「ふんふんなるほどね」と理解できるタイプの人間ではないので、これから1ページずつガシガシ鉛筆で書き込みしながら、丁寧に読んでいこうと思います。例えば110ページを開いてみると、文字の半分くらいが「X」なんですが、もしかしたら僕以外にも、これをみて「ｳｯ...」と感じてしまう人はいるかもしれません。

でも恐らく、本書ただ一冊を読んで完結するようなことは想定されていなくて、「データサイエンス入門シリーズ」の他の本も併せて読んで勉強することが想定されているのだと思います。それらも含めて丁寧に読み進めていけば、きっと怖くないと思うので、頑張って読んでいきたいと思います。

まとめ

決して阿っているわけではなくて、本当にいい教科書だなと思いました。時間をかけて、Rによるプログラミングと数理の勉強を並行していくタイプの書籍だと思いますので、ゼミなどで教科書として利用することが特に相性が良いかもしれません。

最後に...これは完全に個人的な話ですが、本書の査読者（原稿をチェックされた方々？）のお名前を眺めていたら、学部の同期がいました。何よりそれが一番の励みになりました。

2020-12-15

複数の.csvファイル（じゃなくてもいいけど）を一気に結合する

R オンライン実験

本記事について

Online Psychological Experiment Advent Calendar 2020 第17日目の記事です。

最初は気軽にポエムでも書こうかと思っていたものの、他の記事がどれも真面目だったので、「だから、オレだってなんかしなくっちゃあな…カッコ悪くて2021年に行けねーぜ…」ということで急きょtipsを書くことにしました。

オンライン実験のデータがどのように保存されるか

もちろん、どんな実験をするかや、どのサービスを利用するかによって、話は変わります。しかし一般的には、「協力してくださった実験参加者の数だけ、ファイルが蓄積される」と思います。オフラインの（実験室で行うような）実験ならなおさらですね。

例えば僕自身は、こちらの記事を参考にして、lab.jsで作成した.jsonファイルを、Open Labにアップロードし、単語記憶課題のオンライン実験を行いました。その結果、実験参加者の数だけ個別の.csvファイルがサーバ上に蓄積されました。

複数の.csvファイル（じゃなくてもいいけど）を一気に結合したい

もちろん.csvファイルに限らず、全てのファイルが.xlsxでも良いのですが、ともかく拡張子が同じ複数のファイルを、1つのファイルに結合しなければ、前処理も分析もできないわけです。

実験参加者数が少なければ（例えば10人など）、1つのファイルを開いてデータをコピーして、別のファイルにペーストするのを繰り返す...という手作業でも、さほど手間にはならないでしょう（※言うまでもなく、これはヒューマンエラーが混入する恐れがあるので、避けるべきですが）。

しかしオンライン実験となると、一般的にサンプルサイズが大きくなりやすいと思うので、手作業は労力の観点から現実的ではありません（※くどいようですが、労力以外の観点からも、データをコピペでまとめるのは避けるべきです）。

purrrパッケージを使おう

R環境限定の話になりますが、purrrパッケージを用いれば、同じディレクトリ内に存在する、拡張子が同じ全てのファイルを、1つのデータフレームに結合することが、容易にできます。

その方法は...これだっ！（他の方が作成された資料を共有していくツェペリ魂）

上の画像のコードは、tidyverseパッケージをロードする前提で書かれているので、最低限必要なpurrrパッケージとreadrパッケージだけをロードした場合には、以下の書き方になります。
（tidyverseパッケージをロードしたら、これら2つのパッケージも自動的にロードされますが）

library(purrr)
library(readr)


target_files = list.files("data/",
                          pattern = ".csv$",
                          full.names = TRUE)

dat = purrr::map_df(target_files, readr::read_csv)

まず、特定のディレクトリ内に存在する、拡張子が.csvのファイル名を全て、target_filesという名前のリストにまとめています。この場合は、現在の作業ディレクトリの中に、「data」という名前のフォルダがあり、その中に全ファイルが存在する想定です。

次に、 readr::read_csv()関数で個々の.csvファイルを読み込むたび、purrr::map_df()で1つのデータフレームとして結合していきます。結合されたデータフレームは、datという名前のオブジェクトに格納されます。

自分用にカスタマイズも可能

データによっては、

ファイルを読み込む際に、空欄のセルをNAにしたい
ファイル名も、データフレーム内の変数として保存したい

などの、細かいカスタムが必要となることもあるでしょう。そのような場合には、データを読み込む関数（上のコードではreadr::read_csv()）をカスタムした自作関数を定義すればよいです。

ここで新たにdplyrパッケージをロードしているので、パイプ演算子%>%を使用しています（それだったらtidyverseパッケージを1つロードすればよい気もしなくもない）。

library(purrr)
library(readr)
library(dplyr)

# 自作関数
read_func = function(x){readr::read_csv(file = x, na = "") %>%
    dplyr::mutate(filename = x)}


target_files = list.files("data/",
                          pattern = ".csv$",
                          full.names = TRUE)

dat = purrr::map_df(target_files, read_func)

このコードでは、read_func()という関数を自作しています。
まずreadr::read_csv()でファイルを読み込む際に、空白のセルをNAにするよう、引数を指定しています。さらに、1つのファイルをそのようにして読み込んだあとで、dplyr::mutate()によって、ファイル名をfilenameという名前の列に追記しています。
あとはpurrr::map_df()のなかで、自作した関数read_func()を使用することを明記するだけです。

おわり

全てのデータが、1つのデータフレームにまとめられれば、第15日目の記事のように前処理を行うことも、分析することもできます。

Enjoy !!

2020-12-11

Google Colaboratory + PyStanを使ってみた話

Stan Python

この記事について

Stan Advent Calendar 2020 第14日目の記事です。
14日目のカレンダーが空いていたので、何か埋めなきゃということで、「そういえば使ったことないから、Python環境でStanを使ってみよう」という思い付きで書いています。

普段はR環境でStanを使っているので、全てが手探りです。そんなわけで、この記事は自分のための備忘録です。世の中にはもっと体系的にまとまった記事がたくさんあるので、ぜひそれらを参考にしてください。

今回はGoogle Colaboratory上で、PyStanを用いることにします。なお今回はGoogle Driveのマウントはしませんが、Google ColaboratoryをGoogle Driveと連携させることで、より使いやすくなると思います。

実演

.stanファイルを用意する

これはRやPythonとは関係のない作業なので*1、好きなエディタを用いて書きます。ここでは説明変数が1つの、単回帰モデルを例とします。

$y_i=\beta_0+\beta_1x_i+\varepsilon_i \qquad i = 1, 2, ..., n \\ \varepsilon_i \sim {\rm Normal}(0, \sigma)$

$n$ 個のデータがあり、目的変数 $y$ を、説明変数 $x_1$ で予測しています。 $\beta_0$ が切片、 $\beta_1$ が回帰係数、 $\varepsilon$ が残差です。
これは以下の確率モデルと同義です。

$y_i \sim {\rm Normal}(\beta_0+\beta_1x_i, \sigma) \qquad i = 1, 2, ..., n \\$

つまり、

目的変数 $y_i$ は、正規分布に従う。その平均は $\beta_0+\beta_1x_i$ で、標準偏差が $\sigma$ である

というモデルを立てたことになります。以下のStanコードを、regression.stanという名前で保存しておきます。

data {
  int<lower=0> n_data;   //サンプルサイズ
  real y[n_data];        //目的変数
  vector[n_data]  x;     //説明変数
}

parameters {
  real beta_0;           //切片
  real beta_1;           //回帰係数
  real<lower=0> sigma;   //標準偏差
}

model {
  y ~ normal(beta_0 + beta_1*x, sigma);
}

Google Colaboratory上に.stanファイルをアップロードする

自分のPCにPython環境がなくても、インターネットブラウザ上でPythonコードが実行できる、Google Colaboratoryを使用します。ここにアクセスして、[ファイル　→　ノートブックを新規作成]します。

f:id:das_Kino:20201211103021p:plain

ローカルの.stanファイルは、以下の手順で容易にアップロードできます。 f:id:das_Kino:20201211103808p:plain

「注: アップロードしたファイルはランタイムのリサイクル時に削除されます。」というメッセージが現れると思いますが、これはGoogle Colaboratoryを使う以上仕方のないことです。

アップロードされた場所は、対象のファイル名にカーソルを合わせて、右端の「･･･」をクリックすると現れるメニューのうち、「パスをコピー」を選択するとクリップボードに格納されるので、知ることができます。実際にやってみると、/content/regression.stanであることが分かります。

f:id:das_Kino:20201211104119p:plain

ライブラリの読み込み

最低限、これらのライブラリを読み込めば良いと思います。PyStanのサンプリング結果を、トレースプロットやチェインごとの事後分布といった形で可視化するためには、arvizというライブラリが適しているようですが、Google Colaboratoryではまず!pip install arvizでインストールする必要があります。

import numpy as np #サンプルデータを作成する際に、randomモジュールで乱数生成するため
import pystan #stanのPython用インタフェース

# 可視化用ライブラリ ----------------------
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

# PyStanのサンプリング結果の可視化用ライブラリ ------------------
!pip install arviz # arvizというライブラリをまずインストールする
import arviz

サンプルデータの作成

単回帰モデルを真のモデルとして、パラメータの真値を指定してサンプルデータを作成します。

n = 50 #サンプルサイズ

np.random.seed(seed = 1234) #乱数のシード
x = np.random.randint(low = 1, high = 30, size = n) #説明変数。1～29の一様乱数（整数）をn個生成
y = np.random.normal(loc = 20 + 1.5 * x,
                     scale = 3.0,
                     size = n) #目的変数。正規分布に従い、その平均は20 + 1.5x。標準偏差は3.0

可視化してみます。

# サンプルデータの可視化
sns.jointplot(x = x, y = y)

f:id:das_Kino:20201211110044p:plain

.stanファイルのコンパイルとサンプリング

まずはアップロード済みの.stanファイルをコンパイルします。

sm = pystan.StanModel(file = '/content/regression.stan')

次に、サンプリングです。Rだと、rstan::sampling()関数の中で、object = smとコンパイルしたオブジェクトを指定しますが、Pythonだとちょっと書き方が違うことに注意です。

また、Rでは、rstan::sampling()関数の中で、list型でデータを渡しますが、PyStanでは辞書型でデータを渡す必要があるところも違いますね。辞書型のデータはdict()などで作成できます。

その他の引数はRStanでもPyStanでもほとんど同じですね。

fit = sm.sampling(
    data = dict(
        n_data = n,
        y = y,
        x = x
        ),
    seed = 1234,
    iter = 2000,
    warmup = 1000,
    chains = 4)

収束診断とサンプリング結果の出力

まずはうまく収束したかどうか、視覚的に確認してみましょう。ここではarvizライブラリのplot_trace()を用いて、トレースプロットやチェインごとの事後分布を可視化します。

arviz.plot_trace(fit)

f:id:das_Kino:20201211111419p:plain

収束していると考えてよさそうです。
要約統計量は以下の通りです。うまくパラメータリカバリできていますね。

print(fit)

f:id:das_Kino:20201211111202p:plain

おわり

PyStan、初めて使ってみたんですが、Google Colaboratoryを併用することで、環境構築の手間もほとんどいらず、手軽に利用することが出来ました。

Enjoy !!

*1:もっとも、RStudioの文法チェック機能が優れているので、慣れないうちはRStudio上で書くのが良いような気はします

2020-12-09

Stan Advent Boot Camp 第9日目: 階層線形モデルをやってみよう

R Stan

この記事について

Stan Advent Calendar 2020 9日目の記事です。本記事で紹介する階層線形モデルは、第4日目に紹介された重回帰モデルや、第5日目に紹介されたロジスティック回帰モデルなどの、一般線形モデルや一般化線形モデルを拡張したモデルになります。Stanで階層線形モデルの確率モデルを書く方法を紹介します。

なおこの記事では階層線形モデルという名称を使用しますが、線形混合モデルやマルチレベルモデルとも本質的には同じです。

一般線形モデルの復習

ここでは説明変数が1つの、単回帰モデルを例とします。

$y_i=\beta_0+\beta_1x_i+\varepsilon_i \qquad i = 1, 2, ..., n \\ \varepsilon_i \sim {\rm Normal}(0, \sigma)$

$y_i \sim {\rm Normal}(\beta_0+\beta_1x_i, \sigma) \qquad i = 1, 2, ..., n \\$

つまり、

目的変数 $y_i$ は、正規分布に従う。その平均は $\beta_0+\beta_1x_i$ で、標準偏差が $\sigma$ である

というモデルを立てたことになります。

階層線形モデルへの拡張

データの階層性を考慮すべき状況

StanとRでベイズ統計モデリングの121ページから引用すると、階層モデルとは

説明変数だけでは説明がつかない、グループに由来する差異（グループ差）を上手く扱うための一手法

です。
「グループ」というのは例えば、

実験参加者ID	グループ
1	◇◇大学
2	◇◇大学
3	◇◇大学
4	▲▲大学
5	▲▲大学
6	▲▲大学
7	☆☆大学
8	☆☆大学
9	☆☆大学

などです。グループは上表のような「所属集団」に限らず、以下のように個人をグループと見なすことも可能です。下表では、「実験参加者ID」がグループに相当しています。

実験参加者ID	試行番号
1	1
1	2
1	3
2	1
2	2
2	3
3	1
3	2
3	3

上の2つの表ではいずれも、9つのデータが得られていることになっています。では、目的変数 $y_i$ と説明変数 $x_i$ も同様に9つずつデータが得られているとき、

$y_i \sim {\rm Normal}(\beta_0+\beta_1x_{i}, \sigma) \qquad i = 1, 2, ..., n \\$

という回帰モデルを考えてよいかというと...そうとは限りません。ここで、データの階層性を考慮しなければならない場合があります。その理由は、以下のスライドの12枚目に詳しいです。

マルチレベルモデル講習会理論編 from Hiroshi Shimizu

階層線形モデルの確率モデル

上記の回帰モデルを階層線形モデルに拡張したい場合、色々な拡張の方法がありますが、典型的には以下の3通りが考えられます。

切片のみに、グループ差を考慮する
回帰係数のみに、グループ差を考慮する
切片と回帰係数に、グループ差を考慮する

ここでは一例として、「3. 切片と回帰係数に、グループ差を考慮する」場合の階層線形モデルの確率モデルを紹介します。

$y_{i} \sim {\rm Normal}(\beta_{0j\lbrack i \rbrack}+\beta_{1j\lbrack i \rbrack}x_{i}, \sigma) \\ \beta_{0j} \sim {\rm Normal}(\mu_0, \sigma_0) \\ \beta_{1j} \sim {\rm Normal}(\mu_1, \sigma_1)$

添え字 $i$ は一つひとつのデータを、添え字 $j$ は一つひとつのグループを表します。上の表でいえば（一つ目の表）、添え字 $i$ は一人一人の実験参加者を、添え字 $j$ は各大学を表しています。切片パラメータ $\beta_0$ と回帰係数パラメータ $\beta_1$ に添え字 $j$ が付いて、 $\beta_{0j}$ と $\beta_{1j}$ になっているので、グループごとにこれらのパラメータが異なることを仮定しています。

では、グループごとに異なると仮定した切片と回帰係数は、どのように、どれくらいグループごとに異なるのでしょうか。それを反映しているのが、

$\beta_{0j} \sim {\rm Normal}(\mu_0, \sigma_0) \\ \beta_{1j} \sim {\rm Normal}(\mu_1, \sigma_1)$

です。つまりグループごとの切片と回帰係数の違いは、正規分布に従うというモデルを立てたことになります。

実演

まずは、パラメータの真値と真のモデルが分かっている状況からサンプルデータを生成します。

library(ggplot2)
library(rstan)

set.seed(1234)

n_data = 100 #総データ数
n_group = 5 #グループ数
n_each_group = n_data / n_group　#各グループ内の人数 = 20

mu_0 = rnorm(n = n_group, mean = 30, sd = 10) #切片は、平均30、標準偏差10の正規分布に従う
mu_1 = rnorm(n = n_group, mean = 3, sd = 1.5) #回帰係数、平均3、標準偏差1.5の正規分布に従う

x = runif(n = n_data, min = 1, max = 30) #説明変数
y = rnorm(n = n_data, mean = (mu_0 + mu_1*x), sd = 15) #目的変数
group = rep(1:n_group, time = n_each_group) #グループを識別する変数

このようなサンプルデータの作り方は、StanとRでベイズ統計モデリングに載っています。

グループごとに回帰直線を描くと、こんな感じです。確かに切片も傾きも、グループごとに異なっています。

ggplot(data = data.frame(x, y, group), mapping = aes(x = x, y = y, color = factor(group))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

f:id:das_Kino:20201208170634p:plain

この確率モデルを指定してパラメータを推定したら、真値に近くなることが期待できるはずです。やってみましょう。Stanコードは以下です。

方法1：回帰モデルを拡張した書き方

data {
  int<lower=0> n_data;      //データ数
  int<lower=0> n_group;     //グループ数
  
  real y[n_data];           //目的変数
  vector[n_data] x;         //説明変数
  
  int group[n_data];        //グループを識別する変数
}

parameters {
  vector[n_group] beta[2];  //beta[1]が各グループの切片
                            //beta[2]が各グループの回帰係数
  
  real mu[2];               //mu[1]が各グループの切片の平均
                            //mu[2]が各グループの回帰係数の平均
                            
  real<lower=0> sigma[3];   //sigma[1]が各グループの切片の標準偏差
                            //sigma[2]が各グループの回帰係数の標準偏差
                            //sigma[3]が目的変数が従う正規分布の標準偏差
}

model {
  beta[1] ~ normal(mu[1], sigma[1]);
  beta[2] ~ normal(mu[2], sigma[2]);
  
  y ~ normal(beta[1, group] + beta[2, group].*x, sigma[3]);
}

.*という演算子が用いられていることに注意してください。これは同じ長さのベクトルを、要素ごとに掛算するための演算子です。

これをmixed_1.stanという名前で保存して、以下のRコードからサンプリングを行います。

sm_1 = rstan::stan_model("mixed_1.stan")

fit_1 = rstan::sampling(object = sm_1,
                        data = list(
                          n_data = n_data,
                          n_group = n_group,
                          y = y,
                          x = x,
                          group = group
                        ),
                        seed = 1234,
                        iter = 4000,
                        warmup = 1000)

print(fit_1)

トレースプロットやchainごとの事後分布を視覚的にチェックした限りでは、収束していると判断して良さそうです。

f:id:das_Kino:20201209102951p:plain

事後分布の要約統計量は以下の通りです。

f:id:das_Kino:20201209102236p:plain

グループごとの切片が従う正規分布の平均 $\mu[1$ ]は、真値が30のところ、事後中央値は28.27
グループごとの切片が従う正規分布の標準偏差 $\sigma[1$ ]は、真値が10のところ、事後中央値は18.26
グループごとの回帰係数が従う正規分布の平均 $\mu[2$ ]は、真値が3のところ、事後中央値は2.86
グループごとの回帰係数が従う正規分布の標準偏差 $\sigma[2$ ]は、真値が1.5のところ、事後中央値は1.62
目的変数が従う正規分布の標準偏差 $\sigma[3$ ]は、真値が15のところ、事後中央値は15.22

グループごとの切片が従う正規分布の標準偏差 $\sigma[1$ ]だけちょっと真値と事後中央値が離れていますが、おおむね良くパラメータリカバリできているように思います。

方法2：多変量正規分布を用いた書き方

上記の階層モデルは、多変量正規分布を用いて、以下のように書くこともできます。多変量正規分布は、Stan Advent Calendar 2020 第6日目の記事でも紹介されています。

※前述の確率モデルに対して、違うStanコードの書き方ができるということではなく、確率モデルも多変量正規分布を用いた別のものに変わっていることに注意してください

data {
  int<lower=0> n_data;             //データ数
  int<lower=0> n_group;            //グループ数
  int<lower=0> n_slope;            //グループ差を考慮する回帰係数の数（切片含む）
  
  real y[n_data];                  //目的変数
  matrix[n_data, n_slope] X;       //説明変数（切片含む）の行列
  
  int group[n_data];               //グループを識別する変数
}

parameters {
  vector[n_slope] beta[n_group];   //グループごとの、回帰係数（切片含む）のベクトル
  
  real<lower=0> sigma;             //目的変数が従う正規分布の標準偏差
  vector[n_slope] mu;              //各回帰係数（切片含む）の平均ベクトル

  cov_matrix[n_slope] Cov;         //多変量正規分布の分散共分散行列
}

model {
  for(n in 1:n_data){
    y[n] ~ normal(X[n]*beta[group[n]], sigma);
  }
  
  beta ~ multi_normal(mu, Cov);
}

modelブロックのy[n] ~ normal(X[n]*beta[group[n]], sigma)は、行列の掛算によって、目的変数 $y$ が従う正規分布の平均を計算しています。重回帰モデルにおける同様の書き方は、Stan Advent Calendar2020の4日目の記事でも紹介されています。

ポイントは、各回帰係数（切片を含む）がvector[n_slope] beta[n_group]で宣言されているところです。回帰係数の数（切片を含む）の要素を持つbetaというパラメータが、グループの数（n_group）だけあると宣言しています。

modelブロックのbeta ~ multi_normal(mu, Cov)は、各回帰係数（切片を含む）betaが、多変量正規分布に従うことを仮定しています。多変量正規分布のパラメータは、以下の2つです。

平均ベクトルmu
分散共分散行列Cov

要するに、それぞれの回帰係数は同時に正規分布に従っており、かつ互いに相関があるというモデルになっています。

ただし多変量正規分布の分散共分散行列を、quad_form_diag()というStanの関数を用いて以下のように定義すると、サンプリング効率がよくなるようです。

data {
  int<lower=0> n_data;             //データ数
  int<lower=0> n_group;            //グループ数
  int<lower=0> n_slope;            //グループ差を考慮する回帰係数の数（切片含む）
  
  real y[n_data];                  //目的変数
  matrix[n_data, n_slope] X;       //説明変数（切片含む）の行列
  
  int group[n_data];               //グループを識別する変数
}

parameters {
  vector[n_slope] beta[n_group];   //グループごとの、回帰係数（切片含む）のベクトル
  
  real<lower=0> sigma;             //目的変数が従う正規分布の標準偏差
  vector[n_slope] mu;              //各回帰係数（切片含む）の平均ベクトル
  
  
  corr_matrix[n_slope] Omega;      //それぞれの正規分布同士の相関行列
  vector<lower=0>[n_slope] tau;    //それぞれの正規分布の標準偏差
}

model {
  for(n in 1:n_data){
    y[n] ~ normal(X[n]*beta[group[n]], sigma);
  }
  
  beta ~ multi_normal(mu, quad_form_diag(Omega, tau));
  
  Omega ~ lkj_corr(2);             //相関行列の弱情報事前分布。
}

Stan モデリング言語: ユーザーガイド・リファレンスマニュアル（日本語翻訳版）によればquad_form_diag()という関数は、

quad_form_diag(Sigma,tau)がdiag_matrix(tau) * Sigma * diag_matrix(tau)と等価になるように定義されています. ここで, diag_matrix_(tau)は, 対角成分がtauとなり, それ以外が0の行列を返します

という働きを持ちます（上記のStanコードでは、SigmaというパラメータをOmegaという名前で宣言しています）。

これで本当に分散共分散行列が作れることを確かめてみましょう。今回のように、グループ差を仮定している回帰係数が2つの場合（切片と、説明変数の回帰係数の2つ）、以下のような計算になります。

f:id:das_Kino:20201209121459p:plain

対角成分は、標準偏差に相当する $\tau$ の2乗なので、分散を表します。非対角成分は、標準偏差 $\tau$ の積に、相関係数 $r$ を掛けているので、共分散を表します。確かに、分散共分散行列が作れていますね。

それでは上記のStanコードをmixed_2.stanという名前で保存し、以下のRコードを実行してみましょう。使用しているサンプルデータは、mixed_1.stanに渡したものと同じです。
model.matrix()は、引数にformula（回帰式）を指定すると、説明変数を行列形式でまとめてくれる関数です。今回は切片も回帰係数とみなしているので、全ての値が1の説明変数を追加する必要があるため、X = model.matrix(y ~ 1 + x)という書き方にしています（ただしX = model.matrix(y ~ x)でも同じです）。

X = model.matrix(y ~ 1 + x)  # = の左側は大文字のX、右側は小文字のxであることに注意

fit_2 = rstan::sampling(object = sm_2,
                        data = list(
                          n_data = n_data,
                          n_group = n_group,
                          n_slope = ncol(X),
                          y = y,
                          X = X, # = の左右ともに、大文字のXであることに注意
                          group = group
                        ),
                        seed = 1234,
                        iter = 4000,
                        warmup = 1000)

パラメータ数が多いので、トレースプロットやチェインごとの事後分布は掲載しませんが、収束していると判断して良さそうです。事後分布の要約統計量は以下の通りです。

f:id:das_Kino:20201209122814p:plain

グループごとの切片が従う正規分布の平均 $\mu[1$ ]は、真値が30のところ、事後中央値は28.15
グループごとの切片が従う正規分布の標準偏差 $\tau[1$ ]は、真値が10のところ、事後中央値は19.41
グループごとの回帰係数が従う正規分布の平均 $\mu[2$ ]は、真値が3のところ、事後中央値は2.85
グループごとの回帰係数が従う正規分布の標準偏差 $\tau[2$ ]は、真値が1.5のところ、事後中央値は1.67
目的変数が従う正規分布の標準偏差 $\sigma$ は、真値が15のところ、事後中央値は15.20

相変わらず、グループごとの切片が従う正規分布の標準偏差 $\tau[1$ ]だけちょっと真値と事後中央値が離れていますが、おおむね良くパラメータリカバリできていますね。

まとめ

~~いかがだったでしょうか~~
パラメータに階層性を仮定したモデルは、ちょっとだけ書き方が難しくなりますが、基本的には回帰モデルの拡張と考えられます。実際のデータ分析では、階層性を考慮したほうがよい状況も多いので、参考になれば幸いです。

Enjoy!

2020-12-04

昔はHarukara.Rというのがあったんだが、膝に矢を受けてしまってな

雑談

この記事について

ベイズ塾 Advent Calendar 2020 5日目の記事です。このブログは技術ブログ・書評が中心なんですが、今回の記事はポエムです。

広島ベイズ塾は、創設メンバーの一人、某氏（ここでは仮にA氏とする）が広島にいらっしゃったときにできました。詳しい話は、もしかしたら最終日（2020/12/25）に紹介されるのかもしれません。

磯野ー、カラオケしようぜー

しばらくして、A氏が関西に異動されました。カラオケに誘われました。

A氏「カラオケの前に寿司でも食うか」
僕「おっ、いいですねー。あ、でもどこも混んでますねー。ググったら、近くに1件だけ寿司屋ありました、もうそこにしますか」
A氏「せやな」

繁華街の片隅にひっそりと佇むその店内には、ショーケースに色んなネタの切り身が陳列されていた。ネタの名前も、価格も表示されていない。名もなき詩魚（うお）。

僕「（ひそひそ、Aさん、これ、何のネタやと思います？）」
A氏「（ひそひそ、わからへん。赤いな）」

ええい、残さずに全部食べてやる oh darlin 君は誰

そんなこんなで僕らの間には、「寿司を食べてからカラオケに行く」という謎のルーティーンが出来あがったのでした。僕らがよく行っていた寿司屋の名前（はる〇〇）と、カラオケを合わせて、Harukaraの誕生です。二人組音楽ユニットじゃありません~~（いやそれHALCALI）~~

Rの勉強会でもする？

時を同じくして広島では、R旋風が巻き起こっていました。その中心はHijiyama.RというRの勉強会（僕らも準レギュラーでしたけど）。

関西でもやる？　どちらからともなくそう切り出したのは、偶然ではなく必然だったのかもしれない。

そんなわけで、Rの勉強会をしてから寿司を食べてカラオケに行く、Harukara.Rが誕生しました。

まあRの勉強会といいつつ、その実態は統計の勉強会でした。Rを使って計算をしているんだから間違ってない。ここで階層線形モデルとかベイズ（入門編）とか色々学んでゆきました。広島から参加者が来ることもあって、しかもみんな寿司が好きで歌が上手いので、いつしか常連メンバーが増えていきました。

いいぞ、Harukara.R、どうしてそんなに大きくなってしまったんだい？　真面目にやってきたからよ、ね？~~（アリさんマークのひっこｓ）~~

アヒル本の読書会でもする？

そんなこんなで2016年秋。いつものように寿司を食べていた僕らは、「次のHarukara.Rでは何の勉強する？」と雑談していました。

聞くところによると、StanとRでベイズ統計モデリングという本（のちに、アヒル本という通称が生まれる）が発売になるらしい。しかも著者は、A氏が愛読しているブログの著者らしい。

じゃあもう、この本の読書会するっきゃないでしょ。というわけで当初はHarukara.Rのいち企画として読書会するつもりだったんですが、せっかくだから詳しい人や興味を持っている人もたくさん巻き込んでやろう、ということで、より大きな組織に企画をお任せすることになり、そこで生まれたのがOsaka.Stanでした。

あれよあれよという間に、色んな方にご参集いただくことができ、ついにはアヒル本の著者にまでゲストトークしていただくことになり、ここで多くの出会いが生まれ、色んなことを学びました。

なんかもう、どこに行っても同じ人いるよね

気づけば、色んな勉強会が乱立していたわけですが、どこに行っても同じ人と会うわけです。そんなこんなでいつしか、広島ベイズ塾に色んな人が集まるようになって、今に至るという。

色んな人とつながり、色んな方向から学ばせていただき、一人では知り得なかったことを知ってゆく。人と人を繋ぐ、キューピットの矢を膝に受けて、走り出していくんだ、遠い地平へーーー

（完）

本記事について

本書の特長

1. 簡単な事例で、豊富な理論を紹介

2. Pythonコード

まとめ

本記事について

本書の特徴

1. 「可視化という作業」自体の解説

2. 丁寧な事例解説

3. 地理データの描画

4. 訳者の技量

まとめ

本記事について

本書の特徴

1. 最小限のR環境

2. 明確な道筋

3. 「まずは実行しよう、数理はそれからだ」という帯

まとめ

本記事について

オンライン実験のデータがどのように保存されるか

複数の.csvファイル（じゃなくてもいいけど）を一気に結合したい

purrrパッケージを使おう

自分用にカスタマイズも可能

おわり

この記事について

実演

.stanファイルを用意する

Google Colaboratory上に.stanファイルをアップロードする

ライブラリの読み込み

サンプルデータの作成

.stanファイルのコンパイルとサンプリング

収束診断とサンプリング結果の出力

おわり

この記事について

一般線形モデルの復習

階層線形モデルへの拡張

データの階層性を考慮すべき状況

階層線形モデルの確率モデル

実演

方法1：回帰モデルを拡張した書き方

方法2：多変量正規分布を用いた書き方

まとめ

この記事について

磯野ー、カラオケしようぜー

Rの勉強会でもする？

アヒル本の読書会でもする？

なんかもう、どこに行っても同じ人いるよね

1.　簡単な事例で、豊富な理論を紹介

2.　Pythonコード

1.　「可視化という作業」自体の解説

2.　丁寧な事例解説

3.　地理データの描画　

4.　訳者の技量

1.　最小限のR環境

2.　明確な道筋

3.　「まずは実行しよう、数理はそれからだ」という帯