階層ベイズモデルで勝敗データからプロ棋士の強さを推定する - StatModeling Memorandum
前の記事のモデルを若干拡張して、勝敗データから将棋のプロ棋士の強さ（skill）を推定しました。 まず勝敗データとレーティングの値ですが、こちらのサイトを参考にさせていただきました。このようなデータを日々更新していくのには多大な努力と忍耐がないとできません。素晴らしいサイトです。 モデルのBUGSコードは以下のようになりました。 model{ for (i in 1:N.member) { skill[i, 1] <- skill.0[i] for (t in 2:N.year){ skill[i, t] <- skill[i, t-1] + beta[i] + r.skill[i, t-1] } } for (g in 1:N.game) { winner.p[g] ~ dnorm(skill[Winner[g], Year[g]], tau.k[Winner[g]]) loser.p[
累積和を使って計算の無駄を省く（変化点検出の例） - StatModeling Memorandum
メーリングリストでStanにおいて累積和を使って変化点検出を高速化する話がありましたのでメモです。 ここではRにはじめから用意されているNileのデータに対して変化点検出します。プロットすると以下です。 ここでは、ある変化点より左の部分では平均mu_l・標準偏差sigmaの正規分布に従い、右の部分では平均mu_r・標準偏差sigmaの正規分布に従うとします。 すると、変化点は離散値をとるパラメータなので、周辺化消去しなくてはいけません。単純にはif_else関数を使った以下の実装になります。 7, 8行目：範囲をおおまかに指定しています。これは実行時にデータを1000で割ってスケーリングするので、この値になっています。 しかし、この実装は各cpにおいて、normal_log(Y[t], mu_l, sigma)とnormal_log(Y[t], mu_r, sigma)を重複して評価してい
トピックモデルシリーズ 4 LDA （Latent Dirichlet Allocation） - StatModeling Memorandum
このシリーズのメインともいうべきLDA（[Blei+ 2003]）を説明します。前回のUMの不満点は、ある文書に1つのトピックだけを割り当てるのが明らかにもったいない場合や厳しい場合があります。そこでLDAでは文書を色々なトピックを混ぜあわせたものと考えましょーというのが大きな進歩です。さてこの記事の表記法は以下になります。前回のUMの場合と同一です。 右2列は定数については数値を、そうでないものについてはR内の変数名を書いています。データは前の記事参照。 グラフィカルモデルは以下になります（左: LDA, 右（参考）: 前回のUM）。 見ると四角のプレートがまで伸びてきただけです。しかしながらこれが曲者でUMからかなりのギャップがあります。以下の吹き出しの順に説明していきます。 ① ここではハイパーパラメータからディリクレ分布に従って『文書の数だけ』が生成されます。このは以下のような『文
RStanで『予測にいかす統計モデリングの基本』の売上データの分析をする - StatModeling Memorandum
12/22(日)にBUGS/Stan勉強会#2がドリコム株式会社にて催されました。そこで2つ発表をしました。そのうちの1つ「『予測にいかす統計モデリングの基本』の売上データの分析をトレースしてみた」に関する詳細＆補足＆苦労話をここで書きたいと思います。RStanというパッケージでRからStanというMCMCサンプリングソフトを使っています。 最初に発表内容のスライドは以下になります。ざっと見るにはこれで十分です。 『予測にいかす統計モデリングの基本』の売上データの分析をトレースしてみた from . . 以降ではスライドごとに簡単に補足していきます。 予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書) 作者:樋口 知之発売日: 2011/04/07メディア: 単行本（ソフトカバー） まずは元となった書籍の紹介です。時系列解析の第一人者による分かりやすく丁寧に
あてはめの原理・あてはめを実装する計算法・モデル - StatModeling Memorandum
先日@ibaibabaibaiさんからこんなツイートがありましたので自分用のメモとして残しておきます。 これは久保本も大い問題ありだけど、「モデル」「あてはめの原理」「あてはめを実装する計算法」を混同するのはいい加減にやめたほうが。ＧＬＭはモデル、ベイズや最尤法はあてはめの原理、ＭＣＭＣは計算法。— baibai (@ibaibabaibai) 2014, 5月 30 あてはめの原理の中でもベイズと最尤法ではレベルが違うけど、その辺はおくとして、せめていまの３つは区別しよう。そうでないと「計算が正確になったらモデルの欠点が出てきて予測力が落ちた」なんていうのは理解不能になってしまう。— baibai (@ibaibabaibai) 2014, 5月 30 勝手に整理して表にしました。 あてはめの原理最尤法ベイズ計算法EMアルゴリズム 最急降下法 共役勾配法 などMCMC 変分ベイズ（VB;
変分法をごまかさずに変分ベイズの説明をする - StatModeling Memorandum
StanでADVIが使えるようになったので、変分ベイズの基礎は抑えておきたいなぁと思って最近学んでいました。自分向けのメモとして残します。 対数周辺尤度・変分下限・KL情報量 目的は事後分布の最もよい近似となるを求めることです。にはあとで因子分解可能という条件を入れます。 イエンセンの不等式を使って、対数周辺尤度を下から評価すると、 を変分下限と呼びます。任意の関数の関数です。対数周辺尤度はevidenceとも呼ばれるため、変分下限はevidence lower bound、略してELBOとも呼ばれます。対数周辺尤度と変分下限の差は、 となります。これはと事後分布のKL情報量（Kullback-Leiblerdivergence）です。対数周辺尤度がにはよらない、データのみから決まる定数であることを考えると、事後分布の最もよい近似となるを求めることは、変分下限を最大化することに等価になりま
『用量反応試験における患者の割り付けの深層強化学習による最適化』というタイトルで統計関連学会連合大会で発表しました - StatModeling Memorandum
ありがたいことに統計関連学会連合大会の招待講演の依頼がありましたので喜んで引き受けました。たくさんの質問ありがとうございました。 発表資料を共有します。一言で言うと、臨床試験において各患者を各用量にどう割り付けるのが良いかを強化学習を用いて求める方法です。性能が良く第2相試験の効率を大きく改善すると思っています。実際の臨床試験でぜひ使ってほしいですし、そのための協力は惜しみません。 用量反応試験における患者の割り付けの深層強化学習による最適化 by @MatsuuraKentaro 元論文はこちらです（open access）。 資料の方は分かりやすさ重視のため、評価シナリオにexponentialモデルが入っていないです。論文の方は欠点を明確にするために入っています。 ソースコードは以下です。2022/9/6にPythonのgymライブラリに互換性のない変更が入りましたので、Rayライブ
西浦先生らによる実効再生産数の統計モデルを解説＆拡張する試み - StatModeling Memorandum
先日の西浦先生のニコ生の発表を聞いていない人はぜひ聞いてください。 モデルとデータを以下のリポジトリでオープンにしていただいたので、モデルについて僕が分かる範囲内で少し解説を加えたいと思います。 github.com 実効再生産数を推定するコードが2種類ありまして、最尤推定（Maximum Likelihood Estimation, MLE）を使ったMLE版（Sungmok Jungさん作成）と 、ベイズ推定版（Andrei Akhmetzhanovさん作成）があります。どちらもコンセプトはほぼ同じで、実装が若干異なります。この記事では、ベイズ推定版（以降、元コードと呼びます）の流れを簡単に説明し、その後でその拡張を試みます。 ベイズ推定版の流れ 大きく分けて「データの集計」「back projection」「実効再生産数の推定」の3つの部分からなります。 データの集計 まずは日付ごとの
COVID-19 日本国内の潜在的な陽性者数を推定する試み - StatModeling Memorandum
日本国内の潜在的な陽性者数を推定することは有益ですが、簡単ではありません。PCR検査がランダムになっていないことが推定を難しくしています。有症状者が検査されやすいというselection biasがあるからです。この記事ではいくつか仮定を置いて潜在的な陽性者数を推定したいと思います。 仮定 全国民のうち潜在的に陽性になっている割合 この割合は年代によらず一定と仮定します。ここでは と書きます（posはpositiveの略）。例えば0.0001なら日本人約1億2千万人中、おおよそ12000人が潜在的に陽性になっている計算です。 なお、国民の年代別人口の値はこのページの令和2年3月報 （令和元年10月確定値，令和2年3月概算値） （PDF：301KB） の「2019年10月1日現在（確定値）」の総人口 男女計の値を使用しました。 陽性者中の有症状者の割合 若年層で無症状が多いなど、年代で異なる
NIMBLEでノンパラベイズを試したメモ - StatModeling Memorandum
NIMBLEというRのライブラリがあります。BUGS言語風の文法でC++にコンパイルされるタイプの確率的プログラミング言語です。実装されている推定のアルゴリズムはここに書いてあります。MCMCの他にも以下のようなアルゴリズムがデフォルトで実装されており、実行速度もかなり速いです。 particle filter（Sequential Monte Carlo, SMC） 空間モデルのCARモデル（Gaussian Markov Random Fieldsと等価） ノンパラベイズのChinese Restaurant Process（CRP）とStick-breaking Process（SBP） ここでは、ノンパラベイズのCRPとSBPを試してみます。参考にしたのは公式ドキュメントです。 データは以下の以前の記事と同じで、モデルも似たようなものです。 statmodeling.hatenab
ガウス過程シリーズ 1 概要 - StatModeling Memorandum
Stanのマニュアルの「Gaussian Processes」の章を実際に実行しましたので記録を残します。結論から言いますと、Stanでやる場合は回帰はよいですがクラス分類に使おうとすると計算が遅いし収束も悪いです。 まずGaussian Process（以下GPと呼ぶ）とは何ぞやということですがgpml（ぐぷむる？）として有名な次の書籍の1章が分かりやすいです。→Gaussian Processes for Machine Learning これを咀嚼して勝手に補完してまとめたものが以下になります。 GPは教師あり学習の一手法です。教師あり学習では有限のトレーニングデータから関数を作ることになります。関数はありとあらゆる入力の値に対して予測値を返すものです。この関数を決めるにあたり、2つのアプローチがあります。1つめは関数をあるクラス（例えば線形だとか）に限定するものです。しかしこれは採
Gaussian Process Latent Variable Models (GPLVM) を使ってみる - StatModeling Memorandum
PCA（主成分分析）のド発展版に相当する、ガウス過程を用いたGPLVMをRからサクッと使うまでの備忘録です。 GPLVMの説明で分かりやすいのは、以下の統計数理研究所のH26年度公開講座「ガウス過程の基礎と応用」の持橋先生と大羽先生の発表資料です。 [1] 統計数理研究所 H26年度公開講座「ガウス過程の基礎と応用」 (web) 元論文は以下です。 [2] M. K. Titsias and N. D. Lawrence (2010) Bayesian Gaussian Process Latent Variable Model. Thirteenth International Conference on Artificial Intelligence and Statistics (AISTATS), JMLR: W&CP 9, pp. 844-851. (pdf) そして、提唱者で
TensorFlowで統計モデリング - StatModeling Memorandum
とある勉強会で「TensorFlowで統計モデリング」というタイトルで講義をしました。聴衆はPythonユーザが多く、データ量が大きい問題が多そうだったので、StanよりもTensorFlowで点推定するスキルを伸ばすとメリットが大きいだろうと思ってこのようなタイトルになりました。 発表資料は以下になります。 TensorFlowで統計モデリング from . . 発表資料の途中に出てくるtf_tutorial.htmlとmodeling.htmlの内容は、以下のipynbをhtmlで出力したものです（見づらかったらプログラム名のところをクリックしてGitHubに移動して見てください）。 ちなみに僕が紹介しているTensorFlowの書き方はEager Executionではなく、Define and Runのやや古い書き方です。あまり気にしていませんけど。 ここではTensorFlowを
NUTSとADVI（自動変分ベイズ）の比較 - StatModeling Memorandum
RStan2.9.0がリリースされました。今まで{rstan}パッケージのsampling関数を使っていたところを、vb関数に変更するだけでサンプリングのアルゴリズムをNUTSからADVI（Automatic Differentiation Variational Inference）に変更することができます。ADVIはユーザーが変分下限の導出や近似分布qを用意をすることなしに、自動的に変分ベイズしてくれます。得られるアウトプットはNUTSとほぼ同様で近似事後分布からの乱数サンプルです。ウリはスピードです。NUTSもADVIもデフォルトのオプションのまま実行して、NUTSと比べて50倍ぐらいスピードが出ることもあります。 NUTSと同様にADVIは効率的な探索のため偏微分を使っているので、離散値をとるパラメータは使えませんが、やはり同様に離散パラメータを消去すれば実行できます。そして、微分
PythonのSymPyで変分ベイズの例題を理解する - StatModeling Memorandum
この記事の続きです。 ここではPRMLの10.1.3項の一変数ガウス分布の例題（WikipediaのVariational_Bayesian_methodsのA basic exampleと同じ）をSymPyで解きます。すなわちデータが に従い*1、とが、 に従うという状況です。ここでデータ（）が得られたとして事後分布を変分ベイズで求めます。 まずはじめに、上記の確率モデルから同時分布を書き下しておきます。 なので、 となります。 この問題は単純なので事後分布は厳密に求まるのですが、ここでは変分ベイズで解きます。すなわち、事後分布をで近似します。さらにと因子分解可能と仮定します。そして、前の記事の最後の2つの式を使って、とが収束するまで繰り返し交互に更新して求めるのでした。以下ではこれをSymPyでやります。 from sympy import * from sympy.stats imp
統計モデリングで癌の5年生存率データから良い病院を探す - StatModeling Memorandum
概要 2017年8月9日に国立がん研究センターは、がん治療拠点の約半数にあたる全国188の病院について、癌患者の5年後の生存率データを初めて公表しました（毎日新聞の記事）。報告書は国立がん研究センターが運営するウェブサイトからダウンロードできます（ここ）。報告書をダウンロードしようとすると注意点を記したポップアップが表示されます。大切な部分を抜粋すると以下です。 本報告書には、施設別の生存率を表示していますが、進行がんの多い少ない、高齢者の多い少ないなど、施設毎に治療している患者さんの構成が異なります。そのため、単純に生存率を比較して、その施設の治療成績の良し悪しを論ずることはできません。 一般に高齢者が多い病院ほど、進行癌（ステージが進んだ癌）が多い病院ほど、その病院の生存率は下がるわけです。それならば、統計モデリングで年齢と進行度（ステージ）の影響を取り除いて（専門的な言葉で言えば「調
『ベイズ統計モデリング ―R,JAGS, Stanによるチュートリアル 原著第2版―』 John Kruschke著、前田和寛・小杉考司監訳 - StatModeling Memorandum
タイトルの本を頂きました。ありがとうございます。僕は原著を少し読んだことがあり、こちらで非常に評判が高い本です。翻訳にもかかわらず原著とほぼ同じ値段で購入できます。 先にJAGSになじみのない方へ説明しておきますと、JAGSはRコアメンバーの一人でもあるMartyn Plummer氏によってC++で開発されたMCMCソフトウェアです。Rから使うのが多数派ですが、PythonからもPyJAGSによって使うことができます。 複雑なモデルでなければStanより収束が早く、離散値をとるパラメータも使えるため、プログラミングがそんなに得意でない人がベイズ統計モデリングをはじめるには一番向いていると思います。最近、再び活発に開発され始めたようで、先日JAGS 4.3.0がリリースされました。 JAGS 4.3.0 is released https://t.co/3jExabWcPI— Martyn
MCMCサンプルを{dplyr}で操る - StatModeling Memorandum
RからStanやJAGSを実行して得られるMCMCサンプルは、一般的に iterationの数×chainの数×パラメータの次元 のようなオブジェクトとなっており、凝った操作をしようとするとかなりややこしいです。 『StanとRでベイズ統計モデリング (Wonderful R)』のなかでは、複雑なデータ加工部分は場合によりけりなので深入りしないで、GitHub上でソースコードを提供しています。そこでは、ユーザが新しく覚えることをなるべく少なくするため、Rの標準的な関数であるapply関数群を使っていろいろ算出しています。しかし、apply関数群は慣れていない人には習得しづらい欠点があります。 一方で、Rのデータ加工パッケージとして、%>%によるパイプ処理・{dplyr}パッケージ・{tidyr}パッケージがここ最近よく使われており、僕も重い腰を上げてやっと使い始めたのですが、これが凄く使い
しょラーさんのブログ記事「StanでAizu Online Judgeの難易度・習熟度を推定したい」の追加解析 - StatModeling Memorandum
背景やデータはしょラーさんの以下のブログ記事を読んでください。 kujira16.hateblo.jp この記事ではAOJ-ICPCで付加された貴重な難易度の情報をフル活用して、問題の真の難易度の推定と、各ユーザの習熟度の推定を行います。 この問題の難しさは「解いていない問題が、スキップして取り組んでいないのか、解こうとしたけど解けなかったのか区別できない」という点にあります。そこで、元記事にもあったように問題をスキップする確率を導入してモデリングする必要があります。 とはいえ、まずはモデルのヒントになりそうなグラフを作成します。 以下では元記事にあわせて、難易度をdifficulty（StanコードではD）、習熟度をperformance（Stanコードではpf）と表現します。 データの分布の確認 difficultyの分布 横軸にdifficulty、縦軸に問題の数をとったヒストグラム
スパースモデルではshrinkage factorの分布を考慮しよう ～馬蹄事前分布（horseshoe prior）の紹介～ - StatModeling Memorandum
ベイズ統計の枠組みにおいて、回帰係数の事前分布に二重指数分布（ラプラス分布）を設定し回帰を実行してMAP推定値を求めると、lassoに対応した結果になります。また、回帰係数にt分布を設定する手法もあります。これらの手法は「shrinkage factorの分布」という観点から見ると見通しがよいです。さらに、その観点から見ると、馬蹄事前分布が魅力的な性質を持っていることが分かります。この記事ではそれらを簡単に説明します。 なお、lassoそのものに関しては触れません。岩波DS5がlassoを中心にスパースモデリングを多角的に捉えた良い書籍になっているので、ぜひそちらを参照してください。 岩波データサイエンス Vol.5 発売日: 2017/02/16メディア: 単行本（ソフトカバー） 参考文献 [1] C. Carvalho et al. (2008). The Horseshoe Esti
リリース、障害情報などのサービスのお知らせ
最新の人気エントリーの配信
j次のブックマーク
k前のブックマーク
lあとで読む
eコメント一覧を開く
oページを開く
