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

StatModeling Memorandum

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

2次元以上のスポット検出を行う統計モデル

1次元の場合の変化点検出は以下の記事で扱いました。 状態空間モデルでシステムノイズに非ガウス分布(1次元の変化点検出) - StatModeling Memorandum 二つの時系列データの間に「差」があるか判断するには - StatModeling Memorandum 変化点検出のポイント…

統計・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 …

GPy(Pythonのガウス過程用ライブラリ)の使い方

概要 GPyを用いて、サンプルパスの生成、ガウス過程回帰、クラス分類、ポアソン回帰、Bayesian GPLVMを実行しました。自分用のメモです。 参考資料 [1] 公式ページ [2] 公式のチュートリアル [3] Gaussian Process Summer Schoolsの資料 理論的背景は上記の[…

Bayesian GPLVMをStanで実装してみた

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

条件付き独立と有向分離を用いた統計モデルの妥当性チェック

この記事では条件付き独立と有向分離を使い、作成したモデルが背景知識と齟齬がないかチェックする方法を紹介します。以下の本の4章と5.3.1項を参考にしています。 Modeling and Reasoning with Bayesian Networks 僕が持っているのは上の一版なのですが、新…

Python(PyStan)で「StanとRでベイズ統計モデリング」の5.1節を実行する

StanのPythonバインディングであるPyStanが公開されて久しいですが、検索してもあんまり情報がヒットしません。ちょっと寂しいと思ったので、インストールやtraceplotの出力なども含めて、以下の本の5.1節「重回帰」の一部を実行してみました(ステマです)…

ゲルマン先生の「役に立つ統計用語集」

この記事はゲルマン先生(Andrew Gelman)の許諾を得て、Handy statistical lexiconを日本語訳したものです。元記事の用語集は現在も更新中です。英語に抵抗がない人はぜひ元記事を読んで下さい。訳語に関しては親しみやすさを重視し、多くの日本人にあまり…

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

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

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

この記事の続編です。一緒にやろうという人がなかなか現れないので、一人でたたき台を作りました。 目的 目的は機能の割にお得な割安mobile PCを探すことです。mobile PCの厳密な定義はないのですが、ここではディスプレイが12型~14型で重さが1kg前後としま…

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

僕はベイズ統計モデリングをはじめる前(5年ほど前)までは主に機械学習をしていました。その頃は平易な成書はあまりなくて、サポートベクターマシンの理論の難しい本を読んだり、Weka本(当時はこれ)を読みながら実装していたことを思い出します。Pythonで…

階層ベイズモデルとWAIC

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

拡散方程式の解と粒子数の時間変化

拡散現象は幅広く観測されます。この記事では拡散方程式の解と粒子数の時間変化のグラフを備忘録として残します。 拡散方程式は以下です。導出方法はWebにあふれているので検索してください。 ここでは経過時間、は位置を表す座標と経過時間で決まる確率分布…

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

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

蟻本シリーズ 3 スライド最小値

今回は以下の問題を考えます。 長さNの数列x[1], x[2], ..., x[N]と数Kが与えられます。y[i] = min{x[i], x[i+1], ..., x[i+K-1]} (i = 1, ..., N-K+1)として定義される数列y[i]を計算しなさい。 この問題は両端キュー(デック, deque)を用いることで、なん…

蟻本シリーズ 2 ランダムウォーク

今回は以下のランダムウォークの問題を考えます。 I×Jの大きさのグリッドがあります。(1,1)からスタートして、1ターンに上下左右4マスのうち移動できる方向にそれぞれ確率p1,p2,p3,p4で移動します。いくつかのマスには石が置いてあり、通行不可能になってい…

蟻本シリーズ 1 ナップサック問題

「プログラミングコンテストチャレンジブック [第2版]」(通称:蟻本)という本がとてもよかったので、これから3回にわたって統計モデルを絡めた感謝の記事を書こうと思います。 プログラミングコンテストチャレンジブック [第2版] ?問題解決のアルゴリズム…

人口ピラミッドのAge Heapingを階層ベイズで補正する

1週間ぐらい前に以下のツイートがバズっていました。togetterのまとめはこちら。 インドネシアの人口ピラミッド、どうしてこうなったのか自分の年齢を気にしない文化なのか pic.twitter.com/yPcvUCkpD2— やなせ (@ynsitx) 2016年6月16日 このグラフのソース…

「はじめての統計データ分析」 豊田秀樹のメモ

あとがきと6章のあとにあるQ&Aの節が熱い思いに満ちていてオススメです。2.7節「論文・レポートでの報告文例」もユニークです。学生思いの教育者としての一面を垣間見た気がします。 あとがきに書いてあるように、たしかに初級向けの授業で伝統的な統計学と…

Michael Betancourt's Stan Lectureを開催しました

ドワンゴさんに会場提供していただき、2016/6/4にMichael Betancourt's Stan Lectureを開催しました。実はStanの勉強会というのはこれがはじめてではなく、約2年ほど前に催されていたBUGS/Stan勉強会がもととなっています。 また今回はニコ生で放送したので…

Report of Michael Betancourt's Stan Lecture

We are really happy to hold Michael Betancourt's Stan Lecture on June 4 at DWANGO. To tell the truth, this is not the first Stan meeting in Japan. Three BUGS/Stan meetings were held about 2 years ago. DWANGO is a very famous company for it…

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

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

Stanのマニュアル日本語訳プロジェクト

現在、有志でStan 2.9.0のマニュアルを翻訳しています。 https://github.com/stan-ja/stan-ja まったり進行ですが、すでに日本語訳も95ページほどとなっており、言語仕様や関数のリファレンスなどを除けば5~6割ぐらいに達していると思います。 翻訳の稼ぎ頭…

累積和を使って計算の無駄を省く(変化点検出の例)

メーリングリストでStanにおいて累積和を使って変化点検出を高速化する話がありましたのでメモです。 ここではRにはじめから用意されているNileのデータに対して変化点検出します。プロットすると以下です。 ここでは、ある変化点より左の部分では平均mu_l・…

Tweedie分布のパラメータを推定する

@dichikaさんのブログ記事でTweedie分布の存在を知りました。Stanのメーリングリストでも「推定できないの?」という質問は過去にありましたが、多忙のBobさんからは「summing out(離散値をとるパラメータの和をとって消去)すればできるかもねー」という素…

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

R

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

二つの時系列データの間に「差」があるか判断するには

詳しい経緯はこのまとめを参照してください。時間軸でぶった切って各時点で検定を使う手法は、百歩譲って「差があるかどうか」は判定できるかもしれないけど、「どれほど異なるのか」については何も言えない。「どの時刻から異なるか」についても言えるか分…

NUTSとADVI(自動変分ベイズ)の比較

RStan2.9.0がリリースされました。今まで{rstan}パッケージのsampling関数を使っていたところを、vb関数に変更するだけでサンプリングのアルゴリズムをNUTSからADVI(Automatic Differentiation Variational Inference)に変更することができます。ADVIはユ…

変分法をごまかさずに変分ベイズの説明をする

StanでADVIが使えるようになったので、変分ベイズの基礎は抑えておきたいなぁと思って最近学んでいました。自分向けのメモとして残します。 対数周辺尤度・変分下限・KL情報量 目的は事後分布の最もよい近似となるを求めることです。にはあとで因子分解可能…

岩波データサイエンスvol1のいくつかの例題をStanでやってみた

岩波データサイエンスは従来の書籍とは異なり、以下のサポートページの異様な充実がウリの一つです。 https://sites.google.com/site/iwanamidatascience/ また、StanとRStanのはじめの一歩と言える使い方が載っています。また、作図を含めたソースコードがG…

傾向スコアを使ったベイズモデル

傾向スコアについては以下を参考にして下さい。 [1]統計的因果推論(傾向スコア)の勉強会資料をアプしてみた(web) [2]傾向スコアを用いた共変量調整による因果効果の推定と臨床医学・疫学・薬学・公衆衛生分野での応用について(pdf file) [3]傾向スコア…

「使える大学・使えない大学」の事例から考えるアンケートの解析方法

少し前に週刊ダイヤモンドの記事「使える大学・使えない大学」の結果がインターネット上で話題になっていました。具体的には以下のデータです。 引用元はこちら(参考: Googleブックスの書籍を引用するには) 画像の下の方の注意書きにも注目。有効回答数は…

分布から見た線形モデル・GLM・GLMM

久保さんのみどりぼん勉強会もせっかく催されていることだし、それにちなんだ記事を書きたいと思っていました。ここまでいい加減にGLMとGLMMをすっ飛ばして紹介して、さっさとBUGS/Stanのラビリンスパラダイスへいざないたいなぁという心境をスライドにしま…

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

R

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

ガウス過程シリーズ 3 クラス分類(PRML下 Fig 6.12)

今回はGaussian Processで2値クラス分類を行います。2値なのでlogistic linkをかませばOKです。しかしながら、高速化ができなくなります。Stan manualの中にも登場しますがinfer.netの例題の中の「Short Examples: Gaussian Process classifier」にも登場し…

ガウス過程シリーズ 2 高速化&フルベイズ

前回の記事のスピードアップをします。 まずは分散共分散行列をコレスキー分解して multi_normal() から multi_normal_cholesky() を使うようにする方法です。このテの高速化の基本とのことです。コレスキー分解をするメリットはzがi.i.d.から生成される、す…

ガウス過程シリーズ 1 概要

Stanのマニュアルの「Gaussian Processes」の章を実際に実行しましたので記録を残します。結論から言いますと、Stanでやる場合は回帰はよいですがクラス分類に使おうとすると計算が遅いし収束も悪いです。 まずGaussian Process(以下GPと呼ぶ)とは何ぞやと…

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

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

「ベイズ統計の理論と方法」渡辺澄夫のメモ

ベイズ推測を使う人はもちろんのこと、嫌う人にもぜひ一読をすすめたい書籍です。ただし、メインの定理の証明の部分は、代数幾何学の特異点解消定理を使いますし、その他にも複素関数論・経験過程といった知識を要求されます。これらの事前知識に詳しくない…

「基礎からのベイズ統計学―ハミルトニアンモンテカルロ法による実践的入門」豊田秀樹ら

ハミルトニアンモンテカルロ法を以下HMCと略します。この本の内容は主に次の4つに大別されます。 (1~3章・付録A)確率・ベイズの定理・ベイズ推定のおさらい & 伝統的な統計学とベイズ統計学の比較 (4~5章・付録B)事後分布からのサンプリングの詳細(…

「コマンドラインではじめるデータサイエンス」Jeroen Janssensら

この本ではMasonとWiggins(2010)のデータサイエンスの定義に従って解析をすすめていきます。すなわち、(1)データの獲得、(2)データのクレンジング、(3)データの精査、(4)データのモデリング、(5)データの解釈の5ステップです。(5)はコンピュータの出番が少な…

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

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

RStanで『予測にいかす統計モデリングの基本』の売上データの分析をする

12/22(日)にBUGS/Stan勉強会#2がドリコム株式会社にて催されました。そこで2つ発表をしました。そのうちの1つ「『予測にいかす統計モデリングの基本』の売上データの分析をトレースしてみた」に関する詳細&補足&苦労話をここで書きたいと思います。RStanと…

不等間隔の状態空間モデル

日付単位とかでデータを取ることが多いこのご時世、等間隔の状態空間モデルを使うことが多いと思います。しかし、ふと不等間隔の状態空間モデルってどうやるんだろーとつぶやいたところ、ご指導いただきました。いつも大変感謝です。 .@berobero11 細かく等…

状態空間モデルでシステムノイズに非ガウス分布(1次元の変化点検出)

北川本「時系列解析入門」14.4節の例題をやってみましたので途中経過の作業ログとして残します。 まずは失敗例から。 data { int N; real Y[N]; } parameters { real mu[N]; real<lower=0> s_y; real<lower=0> s_mu; } model { for (n in 1:N) Y[n] ~ normal(mu[n], s_y); for </lower=0></lower=0>…

陽に解ける常微分方程式を使ったモデル

今回はデータの背後に簡単な(陽に解ける)常微分方程式で記述できるダイナミクスがあると仮定して、Stanでパラメータの推定を行いたいと思います。 状況として定期的に(例えば一年ごとに)サービスをリリースした場合を考えます。それらのサービスを使う総…

kivantiumさんのブログ記事「アニメキャラのバストサイズとPixiv R-18タグ率の関係」の追加解析

@kivantiumさんの以下のツイートが面白そうすぎて追加解析してみました。特に2つ目のツイートが重要で、これがないと階層ベイズでやってみようという気は起りませんでした。 調査の結果、アニメキャラのバストサイズとPixivでR-18タグが付く割合の相関係数は…

データ解析で割安賃貸物件を探せ!(山手線沿線編)

@housecat442さんのプレゼンにインスパイアされて、某S社様のサイトからスクレイピングさせていただき家賃予測を行いました。目的は広さ・最寄駅・築年や各種設備の割にお得な割安物件を探すことです。首都圏の賃貸物件を全て扱うのは大変なので、まずは山手…

生存時間分析 - ハザード関数に時間相関の制約を入れる

今回のデータは以下のような、1日ごとに得られる観察打ち切りを含む何らかのイベント発生データです(この記事の最後のRコードで作成しています)。 timecens1111……34134130371431481501501500…… time列はイベント発生の時刻、cens列は打ち切りの場合は0, イ…

JAGSで2次元のマルコフ場モデル

WinBUGS, JAGS, Stanの三姉妹を等しく愛する僕としては何使っても同時分布を自由に使えるんだよ、ということで2次元データのマルコフ場モデルをやりました。 まず1次元データに対するCAR modelは、Itoさんのこの記事とこの記事を参照。Rパッケージの{dlm}や{…