StatModeling Memorandum

StatModeling Memorandum

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

2014-01-01から1年間の記事一覧

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

前の記事のモデルを若干拡張して、勝敗データから将棋のプロ棋士の強さ(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なども含みます。技術ブログは誰かのために書く、自分のメモとして書く、というメリット以上のことがあります。僕が思うことを箇条書きにして書きます。 勉強会で黙っていても話しかけられるようになる 昔は勉…

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

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

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

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

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…