ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2014-03-29 winBUGS 関数パラメタ

[][][][][]winBUGS関数パラメタ一覧編

## 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 ¥sqrt{¥frac{¥tau}{2¥pi}} ¥frac{1}{x} e^{-¥frac{¥tau}{2}(¥log{x}-¥mu)^2}; x ¥ge 0
    • 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)
  • 正規分布 ¥sqrt{¥frac{¥tau}{2¥pi}} e^{-¥frac{¥tau}{2}(x-¥mu)^2};-¥infty < x < ¥infty
    • 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分布 ¥Gamma(¥frac{k+d}{2})}{¥Gamma(¥frac{k}{2})k^{¥frac{d}{2}}¥pi^{¥frac{d}{2}}|T|^{¥frac{1}{2}}¥times (1+¥frac{1}{k}(x-¥mu)^T T (x-¥mu))^{-¥frac{k+d}{2}}; -¥infty < x < ¥infty; k ¥ge 2
    • 行列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分布についてごちゃごちゃしてきたときにはこちらを参照

[][][][][]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(英語版)の記事に揃えた(こちら)