とびこみメイル
R
こんさる．
これが唯一の回答ではぜんぜんないんだけど，
まあ，
一番てぬきして書くとこうなる，
の一例．
# 個体数 N.invidual の個体をグループ数 N.group のグループ
# にランダムに分配する
N.group <- 5
N.individual <- 10
counter <- summary(
as.factor(sample(N.group, N.individual, replace = TRUE)),
maxsum = N.group
)
result <- counter[as.character(1:N.group)]
names(result) <- 1:N.group
result <- sapply(result, function(x) ifelse(is.na(x), 0, x))
print(result)
まあ，
初心者は
summary(as.factor(...))
わざとか思いつかないのかも．
だとすると，
こうか?
N.group <- 5
N.individual <- 10
result <- rep(0, N.group)
names(result) <- 1:N.group
for (i in 1:N.individual) {
g <- sample(N.group, 1)
result[g] <- result[g] + 1
}
print(result)
とはいえ，
いつまでも for (...) { ... }
なのも進歩がないよねえ，
ということで．
さてさて，
ゐんばぐすの
多変量正規分布
dmnorm(mu[], tau[,])
と
Wishart 分布
dwish(R[,], Deg.freedom)
の両関数の挙動をさぐる実験開始．
まずは
dmnorm(mu[], tau[,])
から．
BUGS 言語でこのようにベイズモデルを定義し
model
{
for (i in 1:N) {
Y[i, 1:2] ~ dmnorm(mu[], Tau[,])
}
mu[1:2] ~ dmnorm(Zero[], Tau[,])
}
これを動かす R コードはこうなる．
先日
つくった WinBUGS.R
で定義した関数群のおかげで，
こういう実験が格段にやりやすくなった，
というわけだ
(という効果も期待して作った)．
source("WinBUGS.R")
library(MASS) # for mvrnorm()
N <- 1000
set.data(list(
N = N,
Y = mvrnorm(n = N, mu = c(0, 0), Sigma = diag(1, 2)),
Zero = c(0, 0),
Tau = diag(100, 2)
))
clear.param()
set.param("mu", c(1, 2))
post.bugs <- call.bugs(n.iter = 200, n.burnin = 0, n.thin = 1)
print(post.bugs, digits.summary = 3)
とりあえず，
無相関な二変量正規乱数のサンプルの平均を推定させてみる．
問題なのは決めうち分散共分散行列
Tau = diag(100, 2)
の挙動なんだが
……
実行すると
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
mu[1] 0.017 0.003 1.200e-02 1.500e-02 1.700e-02 0.02 2.400e-02 1.006 600
mu[2] 0.018 0.003 1.100e-02 1.600e-02 1.800e-02 0.02 2.400e-02 1.008 240
となった．
推定値の sd が小さいな
……
ということで，
今度は
Tau = diag(0.01, 2)
で計算させてみると
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
mu[1] 0.012 0.318 -0.572 -0.192 0.011 0.240 0.633 1.002 600
mu[2] -0.015 0.323 -0.675 -0.229 -0.021 0.216 0.575 1.007 250
となった．
つまり WinBUGS の
一変量正規分布 dnorm(mu, tau)
で tau
を分散の逆数を指定するのと同じく，
多変量の場合は分散共分散行列の逆行列を指定しないといけない
(逆数の行列でないことに注意)．
次の実験．
多変量正規分布の分散共分散行列の無情報事前分布として使いたい
Wishart 分布の関数
dwish(R[,], Deg.freedom)
の挙動を調べたい．
実験用コードを下のように変更して，
テストデータを相関のある (相関係数 0.7)
二変量正規乱数としてみる．
source("WinBUGS.R")
library(MASS) # for mvrnorm()
N <- 1000
set.data(list(
N = N,
Y = mvrnorm(n = N, mu = c(0, 0), Sigma = matrix(c(2, 1.4, 1.4, 2), 2, 2)),
Zero = c(0, 0),
Tau = diag(0.01, 2),
R = diag(0.1, 2),
Deg.freedom = 2
))
clear.param()
set.param("mu", c(1, 2))
set.param("tau", diag(1, 2))
set.param("var.cov", NA)
post.bugs <- call.bugs(n.iter = 200, n.burnin = 0, n.thin = 1)
BUGS code のほうはこんなかんぢで．
組み込み関数 inverse(m)
は m の逆行列を生成する．
model
{
for (i in 1:N) {
Y[i, 1:2] ~ dmnorm(mu[], tau[,])
}
mu[1:2] ~ dmnorm(Zero[], Tau[,])
tau[1:2, 1:2] ~ dwish(R[,], Deg.freedom)
var.cov[1:2, 1:2] <- inverse(tau[,])
}
で，
これを計算させると正しく分散共分散行列が推定できてるように見えるよね
……
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
mu[1] 0.050 0.048 -0.048 0.018 0.050 0.083 0.138 0.999 600
mu[2] -0.024 0.044 -0.115 -0.054 -0.023 0.008 0.055 1.001 600
tau[1,1] 0.981 0.043 0.903 0.951 0.977 1.009 1.067 1.010 210
tau[1,2] -0.710 0.038 -0.782 -0.734 -0.710 -0.685 -0.643 1.007 350
tau[2,1] -0.710 0.038 -0.782 -0.734 -0.710 -0.685 -0.643 1.007 350
tau[2,2] 1.014 0.044 0.932 0.984 1.015 1.043 1.098 0.999 600
var.cov[1,1] 2.075 0.091 1.908 2.011 2.070 2.136 2.260 1.001 600
var.cov[1,2] 1.453 0.080 1.287 1.398 1.453 1.507 1.611 1.000 600
var.cov[2,1] 1.453 0.080 1.287 1.398 1.453 1.507 1.611 1.000 600
var.cov[2,2] 2.007 0.091 1.830 1.942 2.008 2.064 2.186 1.000 600
deviance 6397.262 3.490 6393.000 6395.000 6397.000 6399.000 6405.000 1.007 600
つまり関数
dwish(R[,], Deg.freedom)
はウィシャート分布ではなく，
その逆行列を生成する逆 Wishart 分布
(inverse Wishart distribution)
である，
ということが確認できた．
