読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

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

R

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

R

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

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

R

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

統計・R・Stan関連の本、用途別のオススメ10冊

年末年始向けに、比較的読みやすい本を中心にオススメします。 統計学入門 色々読んでみましたが、現在決定版と言えるものは存在しないように思えました。個人的には、シグマと積分の復習、場合の数・数え上げの方法、確率、確率変数、確率密度、度数分布と…

情報量規準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) 参考ページはこちら(金子先生のホームページは…

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

R

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