[Stan]生存時間分析のコードと便利なデータセットについて – かものはしの分析ブログ
都内の事業会社で分析やWebマーケティングの仕事をしています。大学・大学院では経済学を通じて統計解析を行うなどしておりました。企業に勤めてからは、機械学習やテキストマイニング、クローリング技術などに関心を持っています。 Twitterアカウント Mr_Sakaue( SKUE ) GitHub 読書メーター ほしいものリスト はじめに 仕事で生存時間分析を使うことは結構あるのですが、マーケティングの良いデータセットがない印象でブログにしにくいと感じていました。また、Stanでの生存時間分析の事例もあまり把握していません。そこで使えそうなデータセットやStanのコードを探して、そのデータに対して生存時間分析を適用してみたいと思います。 目次 ・生存時間分析とは ・生存時間分析で使えるデータ ・生存時間分析をマーケティングで使う際の用途 ・先行研究 ・生存時間分析で使えるデータセット ・Sta
![[Stan]生存時間分析のコードと便利なデータセットについて – かものはしの分析ブログ](https://cdn-ak-scissors.b.st-hatena.com/image/square/a787dca9147676db06027eb58bb689b8fa4eb62b/height=288;version=1;width=512/http%3A%2F%2Fkamonohashiperry.com%2Fwordpress%2Fwp-content%2Fuploads%2F2018%2F08%2Fanimal_kamonohashi.png)
dlmで状態空間モデル その1: 状態空間モデルおさらい & dlmひとまわり - Shuma Walking on
参考 本 萩原先生の青い本: 基礎からわかる時系列分析 ―Rで実践するカルマンフィルタ・MCMC・粒子フィルタ― ブログ/個人サイト 小野さんのCommandeur & Koopman「状態空間時系列分析入門」をRで再現する dlmだけでなく、KFASの方法もある。 Logic of Blueさんdlmの使い方 dlmだけでなく、状態空間モデルの説明などがとても分かりやすかった。 pdf/論文/パッケージなど dlmのVignettes(簡単な使い方の説明) 1. モデルの記述 線形状態空間モデルは以下のように表される。 \[ y_t = F_t x_t + v_t \\ x_t = G_t x_{t-1} + w_t \\ \] v, wは正規分布に従うノイズ。 これをより確率分布の目線に近づけて書くと以下のようになる。 \[ p(y_t|x_t) = N(F_t x_t, V) \;
Ad-Stock効果を推定しつつ回帰を回したい⑤ - 統計コンサルの議事メモ
背景 しつこいようですが、Marketing Mix Modeling（MMM）の話題です。 先日、こんな面白い論文を見つけました。 GoogleのResearcherによるMMMの論文（彼らはMedia Mix Modelingと呼んでいます）なのですが、ヒルの式を用いて広告のShape効果（Carveture効果）を推定するということをやっています。ここでShape効果・carveture効果とは、メディアの露出量に対する目的変数の反応を示す曲線を指すようで、ヒルの式とは： $$ H(x; K, S) = \frac{1}{1 + (\frac{x}{K})^{-S}} $$ であり、$K > 0$や$S > 0$となるパラメータによってLogやSigmoidの形状を表現することができるようです。 ヒルの式によってxがどのような形状となるか、実際に確認してみましょう。まずはヒルの式を以
統計モデリングで癌の5年生存率データから良い病院を探す - StatModeling Memorandum
概要 2017年8月9日に国立がん研究センターは、がん治療拠点の約半数にあたる全国188の病院について、癌患者の5年後の生存率データを初めて公表しました（毎日新聞の記事）。報告書は国立がん研究センターが運営するウェブサイトからダウンロードできます（ここ）。報告書をダウンロードしようとすると注意点を記したポップアップが表示されます。大切な部分を抜粋すると以下です。 本報告書には、施設別の生存率を表示していますが、進行がんの多い少ない、高齢者の多い少ないなど、施設毎に治療している患者さんの構成が異なります。そのため、単純に生存率を比較して、その施設の治療成績の良し悪しを論ずることはできません。 一般に高齢者が多い病院ほど、進行癌（ステージが進んだ癌）が多い病院ほど、その病院の生存率は下がるわけです。それならば、統計モデリングで年齢と進行度（ステージ）の影響を取り除いて（専門的な言葉で言えば「調
Stanの関数を使ってRを拡張して高速化する - StatModeling Memorandum
（C++に自動で変換される）Stanの関数を使ってRを拡張できる機能が、Stan/RStanの2.16で実装開始されて2.17でほぼ完成しました。Rを高速化するためにC++（とRcpp）はあまり書きたくないけれど、Stanの関数なら書いてもいいよという僕得な機能です。この記事ではその方法を簡単に紹介します。 元にした資料はRStanの開発者であるBenさんがStanCon2018で発表したこちらの資料です。 ここでは例として、以下の2つの関数をRで使えるようにしましょう。 1) 機械学習分野でおなじみのlog_sum_exp関数 引数はN個の正の実数 2) データにemax modelという曲線をあてはめた場合の対数尤度を返す関数 引数はデータ（N個のXとYのペア）とパラメータの値 手順は簡単で以下だけです。 functionsブロックだけ書いたstanファイルを用意する R側でrstan
ctsemパッケージに学ぶ、時系列データの構造方程式モデリングのStanによる実装
はじめに Stan Advent Calender 201720日目の記事です。 構造方程式モデリングを使って, 時系列データを扱うことをD-SEM (Dynamic SEM)と呼びます。 ここでは、D-SEMを使って階層ベイズモデルが実行可能なctsemを紹介します。 ベクトル自己回帰の問題：時間の扱い 心理学では1つの変数の推移を記述, 予測することよりも, ある変数の変化がどのような変数で説明されるかといった形式で, 2つ以上の変数の関係を明らかにすることに関心がある場合が多いです。 時系列データで複数の変数間の関係を扱うモデルは種々ありますが, その代表として, ベクトル自己回帰モデルがあります。ベクトル自己回帰モデルでは, 時間を離散値,等間隔のものとする仮定が置かれていて, この仮定を逸脱した場合には, 推定にバイアスが生じることが知られています。 最近, 心理学の中では, 経
Rで様々な表を書く。<br />パラメータの多いStanの結果も美しい表に。 – MrUnadon – Bayesian Statistical Modelings with R and Rstan
まずはTwitterでこの記事をシェアする author: Unadon (見習い飯炊き兵) 動作環境:Mac OS Sierra 10.12.1; R version３.３.１; rstan 2.10.1 はじめに Rを使う時、私はわざわざ表など使わず、”head(data)”なんかでデータの確認を済ませてしまったりすることが多いです。 でも、他人と共有するとなると、見せ方を考えなければならない。それで、面倒だけどエクセルにコピペして…とやっていました。 また、Rのコンソールに収まりきれないデータを確認したいときなども、何らかの手立てを考える必要があるとおもいます。 今回は、そんなデータの確認に使える”表の出力”についてまとめていきます。 本稿のお品書き package{DT}: dataframeをHTMLの表に一発変換(データが大きい場合有用) package{knitr}: dat
mots quotidiens.
下で書いたような離散確率の時系列の他に, 特に社会科学などで, 観測値が ロジスティック正規分布に従っている場合が多くあると思います。 つまり具体的には, 観測値 y ∈ {1,0} (1:生起, 0:非生起)だとして, y 〜 Bernoulli(σ(x)) = Bernoulli (1 / (1 + exp(-x))) x 〜 N(0,σ2) になっているようなモデル。回帰モデルでは x がさらに wTxと回帰になっている場合を考えますが, 議論は基本的に同じです。 これは多項分布の場合はいわゆる対数線形モデルで, 自然言語処理では通常 gradientを計算してL-BFGSやSGDなどの最適化で解くことが多いと思います。 ただし, 最適化の前提となる共変量xが既知ではなく, 学習途中に決まる 潜在変数だったりすると, 最適化してしまうと最初に変な局所解にトラップされて しまい, 学習
科学者のあり方を変える帰納プログラミング - SmartNews Engineering Blog
こんにちは。スマートニュースの高橋力矢です。前回のブログでデータ分析+ゲーム理論を題材として、帰納と演繹をまとめる利点をお伝えしました。なんらかの入力 (e.g., ゲーム理論における利得表) があり、特定のアルゴリズム (e.g., 各プレイヤーの戦略的意思決定) を記述することで出力 (e.g., ナッシュ均衡) を得るアプローチは、ほとんどのソフトウェア・エンジニアが慣れ親しんでいるプログラミングそのものです。つまり多くのエンジニアが手がけるプログラミングの実態は演繹的プログラミングです。ではこの対極に位置する帰納プログラミング (Inductive Programming) はどの程度進歩しているでしょうか。 帰納プログラミングの一分野である確率プログラミング (Probabilistic Programming) は統計学や機械学習との関係が密接で、日本でも利用者の多いStanを
Stanによる混合ガンマ分布の推定 - Qiita
data { int N; real values[N]; int cluster_num; } transformed data { real min_val = min(values); real max_val = max(values); real var_val = variance(values); } parameters { positive_ordered[cluster_num] mu; positive_ordered[cluster_num] k; simplex[cluster_num] ratio; } transformed parameters { real<lower=min_val/var_val> rate[cluster_num]; for (j in 1:cluster_num) { rate[j] = k[j] / mu[j]; } } mode
Michael Betancourt's Stan Lectureを開催しました - StatModeling Memorandum
ドワンゴさんに会場提供していただき、2016/6/4にMichael Betancourt's Stan Lectureを開催しました。実はStanの勉強会というのはこれがはじめてではなく、約2年ほど前に催されていたBUGS/Stan勉強会がもととなっています。 また今回はニコ生で放送したのですが、400人以上の視聴者がいたのは驚きました。この内容で、しかも英語なのに。欲を言うともっと会場にも来てほしかったです。当日の進行は少し前倒しになり、そのために閲覧できなかった方には申し訳ありませんでした。動画を公開しますのでお許しください。 以下では内容を簡単に紹介します。 Hiroki ITO "Dealing with latent discrete parameters in Stan" 一人目は頼み込んで北海道から（自腹で）来てくれました伊東宏樹さんです。動画は次の僕の分と合わせて以下にな
StanとRでレプリカ交換MCMC（parallel tempering） を実行する - StatModeling Memorandum
本記事は発展的な話題です。かつて@Med_KUさんのブログ記事「てさぐれ!!RStanもの」で出てきた例題は局所最適値（local minimum）が多くて、Stanで実行する際も初期値をかなりピシッと決めておかないとダメな例題でした。 しかし、モデルが高次元になってくると最初から“それらしい”初期値なんて分かりようがないし、そもそも理論的にはどんな初期値からはじめても長い間iterationをとれば大域最適値に到達してほしいです。しかし、Stanとは言え、現実的な時間内では局所最適値につかまります。溝が少し深いと出てこれません。そんな状況を打破するための発展的なMCMCの手法の1つに「レプリカ交換MCMC（replica exchange MCMC）」というものがあります。パラレルテンパリング（parallel tempering）としても知られています。 例によってここでは詳しく説明し
累積和を使って計算の無駄を省く（変化点検出の例） - 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)を重複して評価してい
[R] [Stan] で ベクトル ARMA (VARMA) を推定 - ill-identified diary
概要書きたいネタは他にもいくつかあるが, 下準備に時間の掛かりそうなものばかりなので, とりあえず簡単なネタで投稿. Stan は最近使い始めたばかりなので非効率なところがあるかも*1. 時系列データのモデリングでよく使われる, 自己回帰移動平均 (ARMA) 過程のベクトル版である ベクトルARMA (VARMA) モデルの簡単な説明をする VARMA モデルを stan で推定する場合のサンプルコードも紹介アンド解説 (RStan を利用) 目次 概要 目次 VARMA モデルとは VAR から VARMA へ VARMA の反転条件 推定方法 ARMA(1,1) 2次元のVARMA(1,1) VARMA(p,q) 参考文献 関連記事 VARMA モデルとはARMA モデルについては, 専門書では Hayashi (2000), 沖本 (2010) などで取り上げられている*2し, よ
![[R] [Stan] で ベクトル ARMA (VARMA) を推定 - ill-identified diary](https://cdn-ak-scissors.b.st-hatena.com/image/square/add2f01d11205be356c7526b1e7bf031858f2596/height=288;version=1;width=512/http%3A%2F%2Fcdn-ak.f.st-hatena.com%2Fimages%2Ffotolife%2Fi%2Fill-identified%2F20160214%2F20160214193912.png)
変分ベイズの自分向けの説明
StanでADVIが使えるようになったので、変分ベイズの基礎は抑えておきたいなぁと思って最近学んでいました。自分向けのメモとして残します。 ●対数周辺尤度・変分下限・KL情報量 目的は事後分布 の最もよい近似となる を求めることです。 にはあとで因子分解可能 という条件を入れます。 イエンセンの不等式を使って、対数周辺尤度 を下から評価すると、 を変分下限と呼びます。任意の関数 の関数です。対数周辺尤度はevidenceとも呼ばれるため、変分下限はevidence lower bound、略してELBOとも呼ばれます。対数周辺尤度と変分下限の差は、 となります。これは と事後分布 のKL情報量（Kullback-Leiblerdivergence）です。対数周辺尤度が にはよらない、データ のみから決まる定数であることを考えると、事後分布の最もよい近似と
rstanの結果も整形して描画する - 盆栽日記
ひとりアドベントカレンダー9日目です。 昨日と同じ話をrstanでもやってみます。 昨日の例は以下のようなロジスティック回帰でした。 res_glm <- glm(data=MASS::birthwt, low~age+smoke, family="binomial") ですがStanのコードを書きたくない。 ということで、glmmstanを使います。 glmmstanパッケージはかのHADの開発者である清水先生が開発した大便利パッケージです。 詳しくは以下をご覧ください。 http://norimune.net/2675 glmmstanパッケージを使うとこうなります。 devtools::install_github("norimune/glmmstan") library("glmmstan") fml <- low~age+smoke res_glmmstan <- glmmsta
Accelerating Asymptotically Exact MCMC for Computationally Intensive Models via Local Approximations
We construct a new framework for accelerating Markov chain Monte Carlo in posterior sampling problems where standard methods are limited by the computational cost of the likelihood, or of numerical models embedded therein. Our approach introduces local approximations of these models into the Metropolis-Hastings kernel, borrowing ideas from deterministic approximation theory, optimization, and expe
rstanでお手軽にGLMMができる関数，glmmstan（）を作りました。 | Sunny side up!
今回はRの話です。 社会心理学会の方法論セミナーでもGLMMをとりあげましたが，階層ベイズの話も久保先生のトークの中にありました。 GLMMでは，変量効果の数が増えると最尤法だと推定が難しくなったりするので，ベイズ推定のほうが向いています。 しかし，GLMMを直接ベイズ推定してくれる商用ソフトもあまりない（あっても機能が部分的）ので，stanなどのフリーのソフトに頼らざるを得ません。しかし，stanは初心者にはなかなかとっつきにくいので，今回はだれでも簡単にGLMMがベイズ推定できる関数を作ってみました。 実は過去に，同様にGLMMを簡単にstanで走らせてくれるglmer2stanというパッケージを紹介したことがあります（こちら）。しかし，glmer2stanはあとで挙げるようにいろいろ使い勝手が悪いところもあり，自分用に使いやすいものを作ろうと思ったのがはじまりです。 追記： Sapp
R で 状態空間モデル: 状態空間時系列分析入門を {rstan} で再現したい - StatsFragments
前の記事でもリンクさせていただいているが、サイト 「状態空間時系列分析入門」をRで再現する では以下のテキストを {dlm}, {KFAS} で再現されており非常にありがたい。これらのパッケージの使い方については リンク先を読めば困らない感じだ。 自分も勉強のために似たことやりたい、、でも同じことやるのもなあ、、と考えた結果 同テキストの内容 {rstan} を使ってやってみた。 補足 Stan には状態空間表現用の関数 gaussian_dlm_obs ( 利用例 ) があるのだが、自分は使ったことがない。7章までのモデルは全て漸化式で表現されているため、それらを Stan のモデルとして記述した。 状態空間時系列分析入門 作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇出版社/メーカー: シーエ
■Non-Gaussian Noise状態空間モデル - eizoo3010の日記
2015-04-19 ■Non-Gaussian Noise状態空間モデル R Finance 導入 状態空間モデルの観測ノイズやシステムノイズは正規分布で定式化してあることが多く、多くのパッケージもそのような形で定式化されている。一方で実データが従う分布がlong-tailを持つことは珍しいことではなく、そのようなデータをガウス型の状態空間モデルで定式化することには懸念がある。非ガウスノイズをガウス型に近似して計算する手法もあるようだが厳密分布で計算した場合に比べて推定精度が悪くなるようなので、できるだけ非ガウスノイズを用いて計算したい。結局何がしたいかというと、金融研究 第１８巻第１号 要約 非ガウス型状態空間表現による確率的ボラティリティモデルの推定に書いてあることがやりたいというお話。 モデル ここでは以下の形をしたSVモデルを考える。いわゆる離散の文脈でのSVモデル。 はリターン
リリース、障害情報などのサービスのお知らせ
最新の人気エントリーの配信
j次のブックマーク
k前のブックマーク
lあとで読む
eコメント一覧を開く
oページを開く
