StatModeling Memorandum

StatModeling Memorandum

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

Stan

ガウス過程シリーズ 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)事後分布からのサンプリングの詳細(…

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, イ…

久保緑本11章のマルコフ場モデル(空間構造のあるベイズモデル)

空間ベイズモデルの一種であるマルコフ場モデルがやっと分かった気がしますので備忘録を書きます。 きっかけとなったのは以下の動画の復習です。0:13:50~0:16:30と1:14:50~1:32:20の部分です。統数研の伊庭先生がマルコフ場モデルについての講義をしている…

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

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 berobero11 前半の散布図行列に関しては別途記事を書きましたのでそち…

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のサンプリングにお任せ モデルの拡張が簡単 デメリット 計算が遅い。文…

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

周期性のある変数・循環する変数を含むモデリングを実践しましたので紹介します。 スライドは埋め込んで、ソースコードのコピペ&解説をメインにします。 使用したデータは以下。自由に使ってください。 元データ: data.txt 起床時刻だけ抜き出したもの: dat…

Stanのリポジトリにある「BUGS Example」で修行する

たまにはBUGSやStanの勉強法について書きます。 まずは久保先生の緑本の例題(ただし11章を除く)をBUGSやStanで実装するのがhello worldに相当します。 次にThe BUGS Bookをはじめから読みつつ気になったBUGSコードを実際に書いてみるのがよいと思います。 …

Amazon EC2でstanの計算をさせるまでのメモ

最近、自宅計算サーバのメンテに疲れたので重い計算はクラウドで流すことに決めました。Amazon EC2で流すstanが思いのほか快適なのでそれまでの道のりをメモとして残します。AWSの中の人である(違 @yamakatuさんに途中でめっちゃサポートしてもらいました。…

SEMの多重指標モデルとMIMICモデル

豊田先生らの「原因をさぐる統計学」の読了記念に記事にします。データ及びRの{lavaan}を使ったスクリプトは豊田先生の研究室のホームページからダウンロードできます。1992年の本なのにわざわざ2013年に{lavaan}でやった例を追加していて敬服します。今回は…