StatModeling Memorandum

StatModeling Memorandum

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

Stan

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

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

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

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

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

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

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

岩波データサイエンスは従来の書籍とは異なり、以下のサポートページの異様な充実がウリの一つです。 https://sites.google.com/site/iwanamidatascience/ 岩波データサイエンス Vol.1発売日: 2015/10/08メディア: 単行本(ソフトカバー) また、StanとRStan…

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

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

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

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

ガウス過程シリーズ 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コードは以下です。 …

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

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

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

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

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

北川本「時系列解析入門」14.4節の例題をやってみましたので途中経過の作業ログとして残します。 時系列解析入門作者:北川 源四郎発売日: 2005/02/24メディア: 単行本 まずは失敗例から。 data { int N; real Y[N]; } parameters { real mu[N]; real<lower=0> s_y; re</lower=0>…

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

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

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

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

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

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

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

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

ノンパラベイズ(ディリクレ過程)の実装

BUGS bookの11章の8.1節のディリクレ過程の写経です。データは以下のサポートページ(11.8.1: Galaxy clustering: Dirichlet process mixture models)でWinBUGS用のodcファイルで配布されています。 WinBUGSをインストールしていない人のために.RDataにした…

Zero-Inflated Poisson分布を使った来店人数などのモデリング

東京R勉強会(#TokyoR)で「100人のための統計解析 - 和食レストラン編」というタイトルで発表してきました。スライドは以下になります。 100人のための統計解析 和食レストラン編 from . . 前半の散布図行列に関しては別途記事を書きましたのでそちらを参照…

Bayesian Lassoで特徴選択

Stanのマニュアルの「Point Estimation」の章を試しましたので記録を残します。 increment_log_prob関数を使って重回帰をやります。その後、2通りのLassoで特徴選択をします。Stanでやる場合、ロジスティック回帰などにも簡単に組み込めますので拡張性が高い…

SIRモデルからはじめる微分方程式と離散時間確率過程(後編)

前の記事の続きです。 今回はSIRモデルを人数のまま扱い、確率過程で扱います。このことでモデルはより正確になって定量的になりますが、「時間がたったらどうなるのか?」などの定性的な理解は難しくなります。時間に関しては一日ごとに感染者数が発表され…

SIRモデルからはじめる微分方程式と離散時間確率過程(前編)

今年はデング熱やエボラで騒がれました。そのような感染症の伝播によって感染人数がどのように変化するかを表すモデルはいくつかありますが、最もシンプルなものはSIRモデルというものです。Wikipediaの記事はこちら。 総人口をNとして、Sが感受性人口(まだ…

トピックモデルシリーズ 7 DTM (Dynamic Topic Model) の一種

最後はおまけでLDAに時系列を組み合わせた実装を試してみたので紹介します。 今まで「文書」と呼んできたものを「ユーザー」、「単語」と呼んできたものを「アクセスしたWebページ(≒アクション)」と考えます。ユーザーが1日目~31日目までV種類のWebページ…

トピックモデルシリーズ 6 GaP (Gamma-Poisson Model)

この記事の表記は以下です。Wがbag-of-wordsの行列を表すことに注意してください。 右2列は定数については数値を、そうでないものについてはR内の変数名を書いています。データは前の記事参照。 LDAの記事で、『別の視点から見ると、LDAがやっていることは、…

トピックモデルシリーズ 5 PAM (Pachinko Allocation Model)

LDAの不満点の一つとしましては、トピック間の関係性を全て無視しているところです。例えば、「政治」と「経済」なんかは相関ありそうですよね。そういうトピック間の相関を考慮したモデルとしてはCTM(Correlated Topic Model)があります。実はStanのマニ…

トピックモデルシリーズ 4 LDA (Latent Dirichlet Allocation)

このシリーズのメインともいうべきLDA([Blei+ 2003])を説明します。前回のUMの不満点は、ある文書に1つのトピックだけを割り当てるのが明らかにもったいない場合や厳しい場合があります。そこでLDAでは文書を色々なトピックを混ぜあわせたものと考えましょ…

トピックモデルシリーズ 3 UM (Unigram Mixtures)

次にUMを説明します。この記事の表記法は以下になります。 右2列は定数については数値を、そうでないものについてはR内の変数名を書いています。データは前の記事参照。 グラフィカルモデルは以下になります(左: UM, 右: 前回のNB)。 見比べてもらうと分か…

トピックモデルシリーズ 2 NB (Naive Bayes)

このシリーズははじめの2ステップ(NB→UM→LDA)がとっつきにくいですがそこまで理解すれば後のモデルの拡張はそんなに難しくは感じませんでした。そのためNBから順にしっかり理解することが重要と思います。またNBとUMは文書のトピックが与えられているかそ…

トピックモデルシリーズ 1 概要

Stanでトピックモデルを実装するメリット・デメリットについて簡単に触れたいと思います。 メリット 実装がラク。LDAでも30行ぐらい ややこしい推論部分は一切実装しなくてOK。全部Stanのサンプリングにお任せ モデルの拡張が簡単 デメリット 計算が遅い。文…

循環する変数の統計モデリング

周期性のある変数・循環する変数を含むモデリングを実践しましたので紹介します。 スライドは埋め込んで、ソースコードのコピペ&解説をメインにします。 とある病んだ院生の体内時計(サーカディアンリズム) from . . 使用したデータは以下。自由に使って…