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節で取り上げた
MCMC (Markov chain Monte Carlo methods)については, 詳しい情報が書籍 [1]や Webに溢れていますので, そちらを参照ください。 MCMCは事後分布からのランダムサンプリングを得るための道具ですが, パラメータ推定にも使うことができます。 ベイズ統計では, 普通の統計のようにパラメータをただひとつの真の値が存在するとは考えずに確率変数とします。 単純にベイズ統計によってパラメータを推定するには, 多変量の事後分布の期待値や周辺分布の計算, 多次元積分が必要になってしまいますが, MCMCはこれを現実的に行うための計算手法に関する1990年代に生まれたイノベーションです。 Stanは, C++で実装された確率的プログラミング言語です。DSLとも呼べるかもしれません。注目はMCMCサンプラーの中でも, Hamiltonian Monte Carloを
Chapter 1. Introduction Chapter 2. Brief Introduction to Bayesian Statistical Modeling Chapter 3. Introduction to the Generalized Linear Model: The Simplest Model for Count Data Chapter 4. Introduction to Random Effects: Conventional Poisson GLMM for Count Data Chapter 5. State-Space Models for Population Counts Chapter 6. Estimation of the Size of a Closed Population from Capture-Recapture Data C
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してきましたが、まさかの大物登場やいろいろなネタがあって非常に楽しかったです
Twitterで質問があったので，書いておきます。 Stanでは，functions{}ブロックでユーザーが関数を定義できます。 普通の関数は，Rと同じように作ればいいのでそれほど難しくはないというか，違和感はないと思いますが，Stanではサンプリングのための関数も定義できて，それでMCMCもできてしまいます。 以下では，自分で正規分布の式を書いて定義したnormal2という確率分布の関数を作って，それでMCMCをする，ということをやってみます。 まず，普通にStanに最初から入ってる正規分布を使って平均と標準偏差を推定するコードを書いておきます。これをnormal1.stanという名前で保存します。 data{ int N; //サンプルサイズ real y[N]; //データ } parameters{ real mu; //平均値 real<lower=0> sigma; //標準偏
Stan 2.5では常微分方程式が扱えるようになったということなので、マニュアルを参考にためしてみました。 生態学でつかわれるロジスティック式を解かせてみました。式は以下のようなもの: Nは個体数、rは内的増加率、Kは環境収容力（最大個体数）です。 Stanのコードは以下のように。 functions { real[] logistic(real t, real[] y, real[] theta, real[] x_r, int[] x_i) { real dNdt[1]; real r; real K; real N; r <- theta[1]; K <- theta[2]; N <- y[1]; dNdt[1] <- r * N * (K - N) / K; return dNdt; } } data { int<lower=1> T; real N0; real t0; real
今年はデング熱やエボラで騒がれました。そのような感染症の伝播によって感染人数がどのように変化するかを表すモデルはいくつかありますが、最もシンプルなものは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;
リリース、障害情報などのサービスのお知らせ
最新の人気エントリーの配信
j次のブックマーク
k前のブックマーク
lあとで読む
eコメント一覧を開く
oページを開く