ゼロがなくなったポアソン分布のパラメーターを切断ポアソン分布で推定する: Taglibro de H
先日のモデルですが、Stanでは切断ポアソン分布をつかえるので、それをつかえばよいのでした。 マニュアル13.7節を参考にしました。 data { int<lower = 0> M; int<lower = 1> X[M]; } parameters { real<lower = 0> lambda; } model { for (m in 1:M) { X[m] ~ poisson(lambda) T[1, ]; } } generated quantities { int N; { real p = poisson_cdf(0, lambda); // 失敗確率 N = M + neg_binomial_rng(M, (1 - p) / p); } } Rコードです。データは前のエントリと同じです。 fit2 <- stan("truncPois.stan", data = list(
[Stan] 隠れマルコフモデルのテスト1: Taglibro de H
CmdStan 2.24に隠れマルコフモデルが くみこまれたということでテストです。 まず、出力確率は既知という設定で、遷移確率を推定します。 説明はマニュアル26章にありますので、それをみながら試行錯誤です。hmm_marginal関数に与える引数は、出力のそれぞれが得られる確率の対数、遷移確率、初期状態の確率となります。 以下コードです。 結果です。 > fit$summary() > fit$summary() # A tibble: 11 x 10 variable mean median sd mad q5 q95 rhat ess_bulk ess_tail 1 lp__ -666. -666. 1.27 1.08 -668. -665. 1.00 1793. 2028. 2 p[1,1] 0.900 0.901 0.0188 0.0185 0.867 0.929 1.00 2
Stan: ガウス過程をためしてみたメモ: Taglibro de H
Stanでガウス過程をためしてみたメモです。 だいたいStan Modeling Language User’s Guide and Reference Manualの18章を参考にしました。また、 StatModeling Memorandumに、Stanによるガウス過程の記事があります。 ガウス過程シリーズ 1 概要 ガウス過程シリーズ 2 高速化＆フルベイズ ガウス過程シリーズ 3 クラス分類（PRML下 Fig 6.12） Stanのコードは、Stan example-modelsのgp-predict.stanを使用しました。 Rコードは以下のとおりです。 library(ggplot2) library(dplyr) library(rstan) options(mc.cores = parallel::detectCores()) rstan_options(auto_wri
[Stan] 2回測定のN-mixture model, bivariate Poisson distribution: Taglibro de H
[Stan] 2回測定のN-mixture model, bivariate Poisson distribution [統計] ある生物の個体数を測定するとします。調査箇所数がRで、調査箇所iの個体数Niは、平均λのポアソン分布にしたがうとします。各調査地でT回の観測をおこないます。検出確率をpとすると、調査箇所iのt回目の観測数はbinominal(Ni, p)にしたがいます。すなわち Ni ~ Poisson(λ) Yi,t ~ Binomial(Ni, p) こうしたデータからλとpを推定します。 Stanでは離散パラメーターがあつかえないので、Nを周辺化してモデリングしていましたが、multivariate Poisson分布でモデリングできるというので、そちらで実装してみました。ただし、ここではT=2とします。 参考にしたところ Dennis et al. (2015) Com
Stanのgaussian_dlm_obs()によるトレンドモデル: Taglibro de H
以前うまくいかなかった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 2.10をためす: Taglibro de H
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
Stan: RStan 2.7.0で並列化: Taglibro de H
Stan 2.7がリリースされました。RStanがCRANにはいって、並列化に対応したとのことなので、ためしてみました。 以前の常微分方程式をつかったモデルのパラメーター推定のRコードをRStan 2.7標準の並列化に対応させてみます。環境はOS X 10.10.4、R 3.2.1、RStan 2.7.0-1です。 ## データ set.seed(123) r <- 0.1 K <- 1000 ts <- seq(0, 30) N <- vector("integer", length(ts)) N[1] <- 100 for (t in seq_along(ts)[-1]) { N[t] <- rpois(1, N[t -1] + r * N[t - 1] * (K - N[t - 1]) / K) } ## モデル library(rstan) rstan_options(auto_w
状態空間モデルをためす: Taglibro de H
状態空間モデルの勉強。やはり自分でコードを書いてみないと。 まずはPetris and Petrone (2011)にあった例から。 data(Nile) ## dlm ## ## Petris G., Petrone, S. (2011) State Space Models in R ## Journal of Statistical Software 41(4) ## http://www.jstatsoft.org/v41/i04 ## library(dlm) buildNile <- function(theta) { dlmModPoly(order = 1, dV = theta[1], dW = theta[2]) } fit <- dlmMLE(Nile, parm = c(100, 2), buildNile, lower = rep(1e-4, 2)) modNil
MCMCの勉強(1): Taglibro de H
今さら感はあるが、MCMC (Markov Chain Monte Carlo; マルコフ連鎖モンテカルロ)を使えるようになろうと、まずは簡単な例から試してみた。 手始めに、正規乱数から生成した標本の平均と標準偏差を推定してみる。 やはりRを使用。MCMCpackパッケージを あらかじめインストールしておいて、呼び出す。MCMCpack中のMCMCmetrop1R()関数を利用して、メトロポリス法によるMCMC推定をおこなう。 library(MCMCpack) 乱数系列を初期化。 set.seed(1) 平均10、標準偏差3の乱数を1000個生成して、xに入れる。 m <- 10 s <- 3 x <- rnorm(1000, m, s) MCMC推定に使用する関数を用意する。betaは要素数2のベクトル。beta[1]が平均、beta[2]が標準偏差で、betaを推定する。関数の返り値
MCMC再訪(1): Taglibro de H
秋の某統計高座のテキスト作成のため、もう一度MCMCを最初からやってみる。 例題1: 正規分布することがわかっている ある母集団から、X = (5.89, 5.69, 4.71, 6.73, 4.90, 3.34, 4.60, 4.00) という標本が得られたとき、その母集団の平均と標準偏差をベイズ推定する。 コード ## ## Sample 1 ## library(R2WinBUGS) ## data X <- c(5.89, 5.69, 4.71, 6.73, 4.90, 3.34, 4.60, 4.00) ## model model <- function() { for (i in 1:N) { X[i] ~ dnorm(mu, tau); } ## priors mu ~ dnorm(0.0, 1.0E-6); tau ~ dgamma(1.0E-3, 1.0E-3); si
1
リリース、障害情報などのサービスのお知らせ
最新の人気エントリーの配信
j次のブックマーク
k前のブックマーク
lあとで読む
eコメント一覧を開く
oページを開く
StanでZero-infrated negative binomialモデル: Taglibro de H
Stanで隠れマルコフモデル (3) Multistate model: Taglibro de H