StatModeling Memorandum

StatModeling Memorandum

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

R

『用量反応試験における患者の割り付けの深層強化学習による最適化』というタイトルで統計関連学会連合大会で発表しました

ありがたいことに統計関連学会連合大会の招待講演の依頼がありましたので喜んで引き受けました。たくさんの質問ありがとうございました。 発表資料を共有します。一言で言うと、臨床試験において各患者を各用量にどう割り付けるのが良いかを強化学習を用いて…

階層ベイズモデルとWAIC (その2) 積分消去する変数が2個ある場合

この記事では積分消去する変数が2個ある場合において、階層モデルのWAICを算出します。シンプソンの公式を使う場合とモンテカルロ積分を使う場合の二通り行います。 その前に、積分消去する変数が1個の場合を復習します。以前の記事のモデル2の(2)の予測を考…

西浦先生らによる実効再生産数の統計モデルを解説&拡張する試み

先日の西浦先生のニコ生の発表を聞いていない人はぜひ聞いてください。 【8割おじさん西浦教授に聞く】新型コロナの実効再生産数のすべて オンライン講演会生中継/主催:日本科学技術ジャーナリスト会議 モデルとデータを以下のリポジトリでオープンにしてい…

COVID-19 日本国内の潜在的な陽性者数を推定する試み

日本国内の潜在的な陽性者数を推定することは有益ですが、簡単ではありません。PCR検査がランダムになっていないことが推定を難しくしています。有症状者が検査されやすいというselection biasがあるからです。この記事ではいくつか仮定を置いて潜在的な陽性…

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

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

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

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

統計・機械学習・R・Pythonで用途別のオススメ書籍

比較的読みやすい本を中心に紹介します。今後は毎年このページを更新します。 微分積分 高校数学をきちんとやっておけばそんなに困ることないような。偏微分とテイラー展開は大学演習のような本でしっかりやっておきましょう。ラグランジュの未定乗数法のよ…

Tokyo.R#70で「統計モデリングで癌の5年生存率データから良い病院を探す」というタイトルで話しました

以下のイベントで話しました。 第70回R勉強会@東京(#TokyoR): ATND 発表資料は以下です。 統計モデリングで癌の5年生存率データから良い病院を探す from Kentaro Matsuura 前にやった解析において、最終的なモデルにいたるまでのプロセスとモデリングのコ…

逆温度1の事後分布のサンプルからWBICを計算する

この記事は以下のツイートを拝見してやってみようと思いました。 #統計 #Baysian もしも「元論文の式(20)をβ₁=1, β₂=1/log nの場合に適用した公式を使ってWBICを計算すると事後分布のサンプルの違いによる分散が大きくなる」とか「直接逆温度1/log nの事後分…

Stanの関数を使ってRを拡張して高速化する

(C++に自動で変換される)Stanの関数を使ってRを拡張できる機能が、Stan/RStanの2.16で実装開始されて2.17でほぼ完成しました。Rを高速化するためにC++(とRcpp)はあまり書きたくないけれど、Stanの関数なら書いてもいいよという僕得な機能です。この記事…

Osaka.Stan#5で「MCMCサンプルの使い方 ~見る・決める・探す・発生させる~」というタイトルで話しました

先日、以下のイベントで話しました。 『StanとRでベイズ統計モデリング』読書会(Osaka.Stan#5) : ATND 発表資料は以下です。 MCMCサンプルの使い方 ~見る・決める・探す・発生させる~ from Kentaro Matsuura 理論的には事後分布や予測分布の使い方というの…

{Rcpp}と{RcppGSL}を活用した数値積分の例

R

GSLはGNU Scientific Libraryの略で, 広い用途の数値計算向けのC言語のライブラリです。20年前から開発されており、まだリリースされ続けています。個人的な印象としては、めちゃくちゃチューニングされてはいないが、長年の開発のおかげで安定しており、マ…

MCMCサンプルを{dplyr}で操る

R

RからStanやJAGSを実行して得られるMCMCサンプルは、一般的に iterationの数×chainの数×パラメータの次元 のようなオブジェクトとなっており、凝った操作をしようとするとかなりややこしいです。 『StanとRでベイズ統計モデリング (Wonderful R)』のなかでは…

情報量規準LOOCVとWAICの比較

この記事はStan Advent Calendar 2016およびR Advent Calendar 2016の12月7日の記事です。StanコードとRコードは記事の最後にあります。 背景は以下です。 [1] Aki Vehtari, Andrew Gelman, Jonah Gabry (2015). Practical Bayesian model evaluation using …

Bayesian GPLVMをStanで実装してみた

この記事の続きです。PRML下の12章に出てくるOil Flowのデータ(データ点1000個×特徴量12個)に対してBayesian GPLVMで2次元(または3次元)の潜在変数空間にマッピングして綺麗に分離されるか見てみます。 まずはPRMLにもあるように普通の主成分分析でやる…

階層ベイズモデルとWAIC

この記事では階層ベイズモデルの場合のWAICとは何か、またその場合のWAICの高速な算出方法について書きます。 背景 以下の2つの資料を参照してください。[1]に二種類の実装が載っています。[2]に明快な理論的補足が載っています。 [1] 階層ベイズとWAIC (清…

underdispersion(過小分散)な場合のポアソン分布の代替

overdispersion(過分散)なポアソン分布は個体差&ポアソン分布で説明するのがシンプルで解釈しやすくて、個人的には好みです。ただ、個体差を考慮するモデルではunderdispersion(過小分散)の場合に対応できません。そのような場合には「ほぼ確定的な値が…

StanとRでレプリカ交換MCMC(parallel tempering) を実行する

本記事は発展的な話題です。かつて@Med_KUさんのブログ記事「てさぐれ!!RStanもの」で出てきた例題は局所最適値(local minimum)が多くて、Stanで実行する際も初期値をかなりピシッと決めておかないとダメな例題でした。 しかし、モデルが高次元になってく…

散布図行列を描くには (corrplot, pairs, GGally)

R

データが与えられた時にはまず可視化をします。そのデータがどのような仕組み(メカニズム)で作られてそうなったかを考えるために必須のプロセスです。しかしながら、どんな可視化がベストかははじめの段階では分からず、とにかくプロットしまくることにな…

Gaussian Process Latent Variable Models (GPLVM) を使ってみる

R

PCA(主成分分析)のド発展版に相当する、ガウス過程を用いたGPLVMをRからサクッと使うまでの備忘録です。 GPLVMの説明で分かりやすいのは、以下の統計数理研究所のH26年度公開講座「ガウス過程の基礎と応用」の持橋先生と大羽先生の発表資料です。 [1] 統計…

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

前回の理論的なまとめを踏まえてStanでWAICとWBICを計算してみます。 今回は例題として混合正規分布から発生させたデータ100個を用いて、2種類のモデルで推定を行い、それぞれに対してWAICとWBICを求めてみます。まずはデータ生成部分のRコードは以下です。 …

「データ分析プロセス」福島真太朗

書籍のタイトルは「データ分析プロセス」とありますが、偉い人を説得してどのようにデータを集めていくかを決めて、KPIをどう設定して~という、いわゆる啓蒙書ではありません。すでに顧客の行動データやPOSデータなどをデータベースに格納しつつあり、そこ…

階層ベイズモデルで勝敗データからプロ棋士の強さを推定する

前の記事のモデルを若干拡張して、勝敗データから将棋のプロ棋士の強さ(skill)を推定しました。 まず勝敗データとレーティングの値ですが、こちらのサイトを参考にさせていただきました。このようなデータを日々更新していくのには多大な努力と忍耐がない…

広津先生による時系列のクラスタリング手法の実装

R

お正月に広津先生のクラスタリング手法をRで実装しました。個体ごとの時系列データをクラスタリングするのに使えます(実際は時系列ではなく一般の2-wayのデータに適用できます)。 個人的に感じる適正なサンプルサイズと時点のサイズはおよそ、10~1000個体…

JAGSの結果からbugsオブジェクトを作成する

R

JAGSを使うとcoda.samples関数で最終的にmcmc.listが得られます。ここからbugsオブジェクトに変換するスクリプトを作りました。 {rjags}, {R2WinBUGS}の2つのライブラリを読んでおく必要があります。 library(rjags) library(R2WinBUGS) mcmc.list2bugs <- f…

infer.netの例題シリーズ4 Bayesian PCA

今回は「Short Examples: Bayesian PCA and Factor Analysis」をやります。みんな大好きPRMLの「12章 連続潜在変数 | 12.2 確率的主成分分析 | 12.2.3 ベイズ的主成分分析」にバッチリ対応します。 infer.netの元ネタと同一の真のパラメータからデータを作っ…

infer.netの例題シリーズ3 混合二変量正規分布のあてはめ

「Tutorial 6: Mixture of Gaussians」をやります。みんな大好き混合正規分布です。 infer.netの例題は分散小さめな2つの正規分布がかなり離れており、点も混じることなくつまらないのでサンプルデータの作成部分を少しいじりました。まず二変量正規分布が2…

Rで谷本距離(Tanimoto Distance)を算出する

R

RでTanimoto Distanceを算出する関数を作りました。 corと同じように使います。粗く言えばcorは相関係数なので「変化が似ている」、distanceは「絶対値が似ている」、Tanimoto Distanceはその間のイメージです。 tanimoto <- function(x) { if(is.data.frame…

RSRubyでRをRubyから使う

RubyからRを使うRSRubyの使い方を説明します。3年ぐらい前から使っているのですがなかなか高速でいいです。 まずはインストール方法から。環境は以下の通りです。 Windows 7 SP1R version 3.0.1 (2013-05-16) 参考ページはこちら(金子先生のホームページは…

アッシェンフェルターのワイン価格予測式をトレースしてみた

1990年頃にアッシェンフェルターさん(Orley Ashenfelter)がビンテージワインの価格予測式(回帰式)を構築しました。精度はかなりよかったものの、当然ワイン評論家からはフルボッコにされました。という話を以前tokyo.RのLTで聞きました。 LTで時間が足り…