StatModeling Memorandum

StatModeling Memorandum

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

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

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

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

モデリングにも役立つ確率分布の性質(再生性と共役事前分布)

自分が分かりやすいように, 応用しやすいように疑似Stanコードで書きました。 再生性 再生性を使うとモデルをシンプルに書けることがあり、推定のスピードアップにつながります。 以下で用いられる2つの確率変数x1, x2は互いに独立とします。 正規分布 x1 ~ …

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

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

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

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

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

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

「The BUGS Book : A Practical Introduction to Bayesian Analysis」 David Lunn et al.

★★★★★の良書です。WinBUGS, OpenBUGSの作者らによる初の教科書です。登場遅すぎですよ。 非常によくまとまっており、久保先生の緑本の次に読むべき本と言えそうです。買いの一択です。 The BUGS Book (Chapman & Hall/CRC Texts in Statistical Science)作者…

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

R

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

確率変数の変数変換とヤコビアンに慣れる

確率変数の変数変換においてヤコビアンを使う場合を簡単にまとめてみました。 ●Q1. 以下は1.00…が出力されます。なぜですか? sum <- 0 max <- 10000 for(i in 1:max){ sum <- sum + var(rnorm(5)) } print(sum/max) 10000回を十分大きな数とみなせば、これ…

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

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

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

R

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

技術ブログを書いてよかったという話

この記事では技術ブログというのはQiitaやRPubsなども含みます。技術ブログは誰かのために書く、自分のメモとして書く、というメリット以上のことがあります。僕が思うことを箇条書きにして書きます。 勉強会で黙っていても話しかけられるようになる 昔は勉…