Menu carp-1.pdf carp-2.pdf carp-3.pdf carp-4.pdf lee-1.pdf lee-2.pdf lee-3.pdf lee-4.pdf bernoulli.stan linear_model.stan roaches.data.R poisson.data.R Hamiltonian Simulations ham-sim-stepsize-just-right.gif ham-sim-stepsize-too-big.gif ham-sim-stepsize-too-small.gif ham-sim.R

1-1. 配列(array)について 配列は任意の型を要素に持つことのできるコンテナです。 したがって、real と int の配列だけでなく、vector や matrix を要素として持つ配列を作ることができます。 二次元配列は、実際は一次元配列を要素に持つ配列です。 real x[M, N] という配列があったとき、x[m, n] は x[m][n] の簡易的な書き方です。 1-2. 行列(matrix)と配列(array)の違い Stan において、matrix と array は明確に異なる役割を持っています。 matrix は数学の行列を表現するためのもので、array は複数の要素を持つためのコンテナです。 したがって、matrix に対しては行列演算(行列積など)を行うことができますが、array に対してはこれを行うことはできません。 また、matrix は線形代数関数(

Stan 2.10がリリースされました（松浦さんのツイート、ChangeLog）。 記法がいろいろかわっているので、Zero-Inflated Poissonモデルを2.10の記法でかいてみました。 // // Zero-inflated Poisson model // data { int<lower=0> N; // Number of observations int<lower=0> Y[N]; // Number of new seedlings vector<lower=0>[N] X; // Proportion of open canopy } parameters { real<lower=0,upper=1> p; // Probability of presence real beta[2]; // Intercept and coefficient } trans

周期性のある変数・循環する変数を含むモデリングを実践しましたので紹介します。 スライドは埋め込んで、ソースコードのコピペ＆解説をメインにします。 とある病んだ院生の体内時計（サーカディアンリズム） from . . 使用したデータは以下。自由に使ってください。 元データ: data.txt 起床時刻だけ抜き出したもの: data_sleepout.txt 起床時刻から起床時刻の時間: data_out2out.txt 起床時刻の累積時刻: data_cum.txt まずはじめに起床時間の分布のパラメータを求めたStanコードは以下になります（model1.stan）。 data { int<lower=0> N; real Sleep_out[N]; } parameters { real<lower=0, upper=2*pi()> mu; real<lower=0, upper=100

概要 前回の(R) Stan の出力加工方法 - ill-identified diaryのおまけ的な形で書いた. BUGS には DIC を計算する機能があるらしいが, rstan にはないので書いてみた. Spiegelhalter et al. (2002) で提案された DIC, デビアンス情報量規準はベイズ統計でモデルの選択に用いられる指標である. 詳しい説明は元論文に任せて, この記事では簡単な説明だけに留める. Gelman et al. (2013) の Ch. 7, takehiko-i-hayashi.hatenablog.com, あるいは 小西 (2008) にも少しだけ言及がある*1. 2017/6/9 追記 不勉強だったのでよく理解していなかったのだが, DIC は計算に事後分布の平均を利用しているため、正則モデルに対してのみ有効である. ベイズ統計に特有の複雑

この記事は、やたらはてブを稼いでしまった前回の記事の続きです。 ASAのプレスリリース及び声明の中には、確かに「p値に依拠しない新たなアプローチの例」として予測値を重視するアプローチ*5、ベイジアンモデリング、決定理論的アプローチ*6およびfalse discovery rate*7といったものを用いるべき、という趣旨のコメントが入っています。とは言え、重回帰分析とか機械学習のような多変量モデリング（なおかつサンプルサイズも大きい）を伴うテーマならともかく、統計学的仮説検定のようなサンプルサイズも小さい（データも少ない）シチュエーションでどうやるんだよ的な疑問を持つ人も多いのではないかと。 そんなわけで、実際にそれっぽい各種検定の数々をStanによるベイジアンモデリングで代替してみたので、この記事ではその結果をつらつら紹介してみようと思います。テーマは前々回のこちらの記事の1節で取り上げた

You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session. You switched accounts on another tab or window. Reload to refresh your session. Dismiss alert

1週間ぐらい前に以下のツイートがバズっていました。togetterのまとめはこちら。 インドネシアの人口ピラミッド、どうしてこうなったのか 自分の年齢を気にしない文化なのか pic.twitter.com/yPcvUCkpD2— やなせ (@ynsitx) 2016年6月16日 このグラフのソースは恐らく以下の西文彦氏による文書です。 統計局ホームページ/統計に関する国際協力/国際協力・交流/インドネシア中央統計庁（BPS）に対する技術協力 Age Heapingの解説がありますので引用します。 Age Heaping（エイジ・ヒーピング）とは、図1に見られるように、人口ピラミッドなどの年齢各歳別の統計において、例えば50 歳、55歳など、0または5で終わる年齢において、人口が突出して多くなる現象のことである。 この現象が起こる原因は、その国において、自分の年齢を正確に知らない人が多いこと

3ヶ月空いてしまいました、が復活します。 あまりに間隔が空いてしまったので、ここ1,2ヶ月の記事を中心にまとめています。ご了承ください。 R3.3.0リリースとR3.3.1プレリリース 上記バージョンが5月上旬にリリースされました。主な変更点についてはこちらをご覧ください。 また、次期バージョンのR3.3.1について、プレリリースバージョンが2016/6/11から、最終リリースバージョンが2016/6/21に予定されています。早すぎる…。 R関連イベント 第53回R勉強会＠東京（#TokyoR）(終了・資料など) 4/30に都内で開催されました。当日の発表で公開されたスライドなどのまとめはこちらでまとまっています: – 第53回R勉強会@東京で発表してきた – INPUTしたらOUTPUT! 私も久々に参加してLTしてきましたが、まさかの大物登場やいろいろなネタがあって非常に楽しかったです

今年はデング熱やエボラで騒がれました。そのような感染症の伝播によって感染人数がどのように変化するかを表すモデルはいくつかありますが、最もシンプルなものはSIRモデルというものです。Wikipediaの記事はこちら。 総人口をNとして、Sが感受性人口（まだ感染してないけど感染する可能性がある人）、Iが発症感染者人口、Rが除外人口（Iから治ってもう感染しない人）を表します。現時点のS,I,Rが与えられた時に、○○時間経った後のS,I,Rはどのようになるでしょうか？○○時間が短い極限をとって連続時間で変化のルールを記述したものを微分方程式と呼びます。SIRモデルでは次のような変化のルールに従うと仮定します。 上記のルールは化学反応式とほとんど同じです。SとIが接触したらある反応速度係数に従ってSもIになる。接触頻度はS*Iに比例するはず。I自体も単体で、ある反応速度係数に従ってRになる。これらを

以前うまくいかなかったStanのgaussian_dlm_obs()によるトレンドモデルだが、モデルを見直してできるようになった。 ナイル川のデータをつかった。Stanによるモデルは以下のように。 data { int<lower=0> N; matrix[1, N] y; } transformed data { matrix[2, 1] F; matrix[2, 2] G; vector[2] m0; cov_matrix[2] C0; F[1, 1] <- 1; F[2, 1] <- 0; G[1, 1] <- 1; G[1, 2] <- 1; G[2, 1] <- 0; G[2, 2] <- 1; m0[1] <- 0; m0[2] <- 0; C0[1, 1] <- 1.0e+7; C0[1, 2] <- 0; C0[2, 1] <- 0; C0[2, 2] <- 1.0e+7;

メーリングリストで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)を重複して評価してい