2014-03-29

winBUGS関数パラメタ一覧編

このマニュアルの57ページからに対応しています
ベルヌーイ分布 $p^r(1-p)^{1-r};r=0,1,...$

## winBUGS
r ~ dbern(p)

#### R
N <- 1000
p <- 0.3
R <- sample(0:1,N,replace=TRUE,prob=c(1-p,p))
hist(R)

二項分布 $\frac{n!}{r!(n-r)!}p^r(1-p)^{n-r}; r = 0,1,..,n$

## winBUGS
r ~ dbin(p,n)

#### R
n <- 100
R <- rbinom(N,n,p)
hist(R)

カテゴリカル $p=(p_1,p_2...,p_n); r = 1,2,..,n; \sum_{i=1}^{n} p_i = 1$

## winBUGS
r ~ dcat(p[])

#### R
n <- 5
p <- 1:n
p <- p/sum(p)
R <- sample(1:n,N,replace=TRUE,prob=p)
tabulate(R,n)

負の二項分布 $\frac{(x+r-1)!}{x!(r-1)!}p^r(1-p)^{x}; x = 0,1,..$

## winBUGS
x ~ dnegbin(p,r)

ポアソン分布 $e^{-\lambda}\frac{\lambda^r}{r!};r=0,1,...$

## winBUGS
x ~ dpois(lambdda)

#### R
lambda <- 3
R <- rpois(N,lambda)
plot(0:max(R),tabulate(R+1),type="h")

ベータ分布 $p^{a-1}(1-p)^{b-1} \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} ; 0 < p < 1$

## winBUGS
p ~ dbeta(a,b)

#### R
a <- 3;b<- 6
R <- rbeta(N,a,b)
hist(R)

カイ二乗分布 $\frac{2^{\frac{k}{2}}x^{\frac{k}{2}-1}}{\Gamma(\frac{k}{2})} e^{-\frac{x}{2}}; x \ge 0$

## winBUGS
x ~ dchisqr(k)

#### R
k <- 3
R <- rchisq(N,k)
hist(R)

Double exponential分布 $\frac{\tau}{2}e^{-\tau |x-\mu|); -\infty < x < \infty$

## winBUGS
x ~ ddexp(mu,tau)

指数分布 $\lambda e^{-\lambda x}; x \ge 0$

## winBUGS
x ~ dexp(lambda)

#### R
lambda <- 2
R <- rexp(N,lambda)
hist(R)

ガンマ分布 $\frac{\mu^r x^{r-1}e^{-\mu x}}{\Gamma(r)}; x \ge 0$

## winBUGS
x ~ dgamma(r,mu)

#### R
r <- 2
mu <- 3
R <- rgamma(N,shape=r,rate=mu)
R. <- rgamma(N,r,mu)
plot(sort(R),sort(R.))
abline(0,1,col=2)

Generalized Gamma $\frac{\beta}{\Gamma(r)}\mu^{\beta r} x^{\beta r-1} e^{-(\mu x)^{\beta}}; x \ge 0$

## winBUGS
x ~ gen.gamma(r,mu,beta)

Log-normal
- winBUGSの $\tau$ とRのsdlog( $s$ )とは $\tau=\frac{1}{s^2}$ の関係であることに注意

## winBUGS
x ~ dlnorm(mu,tau)

#### R
mu <- 3
tau <- 4
R <- rlnorm(N, meanlog = mu, sdlog = 1/sqrt(tau))
hist(R)

Logistic分布 $\frac{\tau e^{\tau(x-\mu)}}{(1+e^{\tau(x-\mu)})^2}; -\infty < x < \infty$

## winBUGS
x ~ dlogis(mu,tau)

正規分布
- winBUGSの $\tau$ とRのsd( $s$ )とは $\tau=\frac{1}{s^2}$ の関係であることに注意

## winBUGS
x ~ dnorm(mu,tau)

#### R
mu <- 2; tau <- 3
R <- rnorm(N,mu,1/sqrt(tau))
hist(R)

パレート分布 $\alpha c^{\alpha} x^{-(\alpha+1)}; x > c$

## winBUGS
x <- dpar(alpha,c)

スチューデントのt分布 $\frac{\Gamma(\frac{k+1}{2})}{\Gamma(\frac{k}{2})}\sqrt{\frac{\tau}{k \pi}} (1+\frac{\tau}{k}(x-\mu)^2)^{-\frac{k+1}{2}}; -\infty < x < \infty; k \ge 2$

## winBUGS
x ~ dt(mu,tau,k)

- Rでは標準正規分布に対応するt分布の関数としてdt(),rt()などがあり、non-central parameterなどを使って対応づけられるが、単純な対応関数にはなっていない
一様分布 $\frac{1}{b-a}; a < x < b$

## winBUGS
x ~ dunif(a,b)

#### R
a <- 2; b <- 5
R <- runif(N,a,b)
hist(R)

ワイブル分布 $v \lambda x^{v-1}e^{-\lambda x^v}; x > 0$

## winBUGS
x ~ dweib(v,lambda)

#### R
v <- 1.5; lambda <- 2
R <- rweibull(N, shape=v, scale = lambda^(-1/v))
hist(R)

多項分布 $\frac{\sum_{i=1}^n x_i)!}{\prod_{i=1}^n x_i!} \prod_{i=1}^n p_i^{x_i}; \sum_{i=1}^n x_i=nx; 0 < p_i < 1; \sum_{i=1}^n p_i=1$

## winBUGS
x[] ~ dmulti(p[],nx)

#### R
n <- 5
p <- runif(n)
p <- p/sum(p)
nx <- 20
R <- rmultinom(N, nx, prob = p)

Dirichlet分布 $\frac{\Gamma(\sum_{i=1}^n \alpha_i)}{\prod_{i=1}^n \Gamma(\alpha_i)} \prod_{i=1}^n p_i^{\alpha_i-1}; 0 < p_i < 1; \sum_{i=1}^n p_i=1$

## winBUGS
p[] ~ ddirch(alpha[])

#### R
library(MCMCpack)
alpha <- c(10,5,3,6)
R <- rdirichlet(N,alpha)

multivariate Normal $(2\pi)^{-\frac{d}{2}} |T|^{\frac{1}{2}} e^{-\frac{1}{2}(x-\mu)^T T (x-\mu)}; -\infty < x < \infty$

## winBUGS muはベクトル、Mは正方行列
x[] ~ dmnorm(mu[],M[,])

#### R
mu <- c(1,2)
sigma <- matrix(c(4,2,2,3), ncol=2)
M <- solve(sigma)
# ただし、winBUGSのMとsigmaの関係はM=solve(sigma)
R <- rmvnorm(N, mean=mu, sigma=solve(M), method="chol")
plot(R)

Multivariate Student-t分布
- 行列Tが分散共分散行列の逆行列になっている点は、multivariate normal distributionの場合と同様
Wishart分布 $|R|^{\frac{k}{2}}|x|^{\frac{k-p-1}{2}}e^{-\frac{1}{2}Tr(Rx)};$ xは対象でpositive definite

## winBUGS
x[,] ~ dwish(R[,],k)

- 多変量正規分布とそのBUGSで使うWishart分布についてごちゃごちゃしてきたときにはこちらを参照

2014-03-29

winBUGS関数パラメタ説明編

BUGS ベイズ MCMC openBUGS winBUGS

昨日の記事でR+openBUGSのことを書いた
winBUGS ~ openBUGSでの分布関数を使ってモデル指定しないといけない
分布関数には形状パラメタとか尺度パラメタとかがあるけれど、場合によって、使い方が違ったりするので、BUGSのモデル設定のパラメタがどういう指定になっているかを知らないと失敗する
- なので、それをきちんと把握するにはBUGSマニュアルで確認するのがよいのだが！
- 色々困ったことがある
  - winBUGSのマニュアルにページ番号がついていないことがある！
  - winBUGSのマニュアルPDFの冒頭目次からクリッカブルになっていない！
- というわけでページ番号のついた版がこちら
- そしてそのp57に関数の引数説明がある→その一覧はこちらの記事
- さて、どうしてそれを確認しないとまずいかと言うと、よく使うdnorm(),dgamma()だけでも要注意だから。
- たとえば、BUGSの正規分布関数 dnorm()では、平均m、分散vの場合には

## winBUGS
dnorm(m,1/v)

- - と書く。これは、Rで値xでの確率密度を

#### R
rnorm(x,m,sqrt(v))

- - と書くのと、紛らわしい
- また、ガンマ関数は正の値だけれど、どんな値をとるか皆目見当がつかないときにBUGSでは

## winBUGS
dgamma(0.0001,0.0001)

- - を使うのがデフォルト(0.0001は『小さい値』のこと)だが、これは、

## winBUGS
dgamma(alpha,beta)

- - として $\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}$ という密度関数に用いられるshapeパラメタとrateパラメタによる指定であり
  - それは、平均 $\frac{\alpha}{\beta}=1$ で、分散 $\frac{\alpha}{\beta^2} =$ とても大、という分布を指定していることになる
  - これのR対応は

#### R
dgamma(x,shape=alpha,rate=beta)

- - Rのdgamma()などでは、紛れの内容にshape,rateを指定して使っておくのが無難（いろんなところで何かしら紛らわしい説明などがあるので）。ここでのshape,rate,scaleパラメタはWiki(英語版）の記事に揃えた（こちら)