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

StatModeling Memorandum

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

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

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

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

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

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

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

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

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

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

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の作者らによる初の教科書です。登場遅すぎですよ。 非常によくまとまっており、久保先生の緑本の次に読むべき本と言えそうです。買いの一択です。 目次はここから。サポートページやサンプルコードはここから。 内容は理…

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

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

ブック40とScanSnap(S1500)で自炊

購入したのはScanSnap S1500 FI-S1500-Aとブック40という裁断機です。 もうネット上に色々な方がレビューしていますが、なおあまり見ない情報というのを参考までにまとめました。以下は処理の順番で書いています。 裁断編 購入動機 ブック40を選んだ動機は価…

独自配列の格安分離型キーボード

ノートパソコンで長時間労働すると肩がこるなど疲労がたまります。そこで分離型キーボードです。Kinesis社のエルゴノミクスキーボードなど魅力的なキーボードはいくつかありますが、以下のようなデメリットがあります。 値段が高い 手がでかい外国人用なので…

所有している理学書、プログラミングの本、推理小説

質問に真摯に答えるために記事を書きます。所有している理学書、プログラミングの本、推理小説を全部載せます。 生物系の本、プログラミングの本はほぼ裁断してpdfにしてしまいました。どうせ全文検索するので、タイトルの入力はあまり頑張ってないです。 推…

PyMC3の計算でGPUを使っている気がしただけの話

たまには浮気させてください。PyMC3は内部でTheanoを使っており、自動微分(auto-diff)が計算可能でStanのサンプラーであるNUTSも実装済みです。またTheanoがGPUに対応しているため、これはMCMCの超高速化が簡単にできるのではッ!と試した記事になります。…

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

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

SEMとベイジアンネットワークとBUGS/Stanの関係

自分用のメモです。 SEM Structural Equation Modelingの略で「構造方程式モデリング」と和訳されています。共分散構造分析は構造方程式モデリング(SEM)の古い名称です。共分散以外の行列も扱えるので、共分散構造分析という名称は適当ではないということ…

あてはめの原理・あてはめを実装する計算法・モデル

先日@ibaibabaibaiさんからこんなツイートがありましたので自分用のメモとして残しておきます。 これは久保本も大い問題ありだけど、「モデル」「あてはめの原理」「あてはめを実装する計算法」を混同するのはいい加減にやめたほうが。GLMはモデル、ベイ…

infer.netの例題シリーズ6 項目反応理論(IRT)+α

今回は「Difficulty versus ability」をBUGSで実装します。この例題のもとはDARE modelというもので以下の論文になります。 "How To Grade a Test Without Knowing the Answers --- A Bayesian Graphical Model for Adaptive Crowdsourcing and Aptitude Tes…

infer.netの例題シリーズ5 Click model (改)

今回は「Short Examples: Click model」をやります。実はこの分析の強化版が「Longer Examples: Click through model」になっていまして、個人的には今回のデータだけでは教育用途かなと思いました。 とにかくデータと目的を紹介します。まずはデータですが…

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…

infer.netの例題シリーズ2 Tutorial 5: Clinical trial

今回は「Tutorial 5: Clinical trial」をやります。データは以下の通り。 outcometreat00001000001101111111 2列目が薬の処理の 有/無、1列目はその結果で 1が治った/0が治らなかった となっています。問題は薬の処理は効果があるかどうか、すなわち治る確率…

infer.netの例題シリーズ1 Bayes Point Machine

BayesPointMachine(以下BPM)の元論文はこちら(R.Herbrich, T.Graepel, and C.Campbell, BayesPointMachines, JMLR, 2001)。 BPMはクラス分類を行うアルゴリズムでSVMに似ています。SVMはマージンを最大にするwを求めますが、BayesPointMachineはモデルの…

WinBUGSのTrapを止めるには

はじめに モデルがよくない/まちがっていることが多いです。徹底的に可視化してからモデル式を考えましょう。 観測データにモデリングを破綻させるような凄い「外れ値」が含まれているとダメ。 全個体の説明変数が同じ値だとダメ。 Linux上でwineを使ってWin…

不等号を含むデータの解析

以下のようなデータがあるとします。 PersonIDYTreat138.3010152.30140.8010226.402020230.102250228.50345.2131357.31359.91467.21448.71460.8147115104.215851591.71561.31572.11 1列目はある人のID、2列目はその人から測定された値(ここでは便宜的にある…

ドラゴンボールの勝敗結果から戦闘力を推定する

以下のような個人データと勝敗データを扱います。 CIDnamepower.level1ベジータ180002悟空80003ナッパ40004クリリン17705栽培マン1200 winnerloser121315233445 モデルを記述する部分であるBUGSコードは以下になります。 model{ for (game in 1:N.game) { wi…

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で時間が足り…

高校生の時に作った傑作詰将棋

高校生の時に作った傑作詰将棋をアップします。僕の棋力はアマ4-5段です。ともに1X手詰めです。腕に自信のある方はチャレンジしてコメントしてくださると嬉しいです。 当時はコンピュータもまともに使えませんでしたが、今はソフトで簡単に余詰チェックとか…