ryamadaのコンピュータ・数学メモ

2013-09-12

ベータ分布を合わせる

いくつかのベータ分布を重み付きで合わせたときの期待値と分散が知りたいとする
単一のベータ分布 $beta(a,b)$ の期待値は $m(a,b)=\frac{a}{a+b}$
単一のベータ分布の分散は $v(a,b)=\frac{ab}{(a+b)^2(a+b+1)}$
$v(a,b)=\int beta(a,b)(x-m(a,b))^2 dx = \int beta(a,b)x^2dx -2m(a,b)\int beta(a,b)x dx +\int m(a,b)^2 dx$ であるが、 $\int beta(a,b) x dx=m(a,b)$ なので
$v(a,b)= \int beta(a,b)x^2dx - m(a,b)^2$
ここに $v(a,b)=\frac{ab}{(a+b)^2(a+b+1)}$ を使えば
$\in beta(a,b)x^2 dx = \frac{a(a+1)}{(a+b)(a+b+1)}$ とわかるから
任意の値の周りの２次モーメントは
$\int beta(a,b)(x-t)^2 dx = \frac{a(a+1)}{(a+b)(a+b+1)}-2t\frac{a}{a+b} +t^2$

# ベータ分布の任意の値を中心とした２次モーメント
my.beta.2moment <- function(a,b,z){
	a*(a+1)/(a+b)/(a+b+1)-2*z*a/(a+b)+z^2
}
# ディリクレ分布に拡張しておく
my.dirichlet.mean <- function(as){
	as/sum(as)
}
my.dirichlet.var <- function(as){
	sum.as <- sum(as)
	as * (sum.as-as)/(sum.as^2 * (sum.as+1))
}
my.dirichlet.2moment <- function(as,t){
	m <- my.dirichlet.mean(as)
	my.dirichlet.var(as)+m^2-2*m*t+t^2
}
#　適当に併せて、乱数を作ってみる
as <- sample(1:10,5)
bs <- sample(1:10,5)
library(MCMCpack)
rr <- rdirichlet(1,rep(1,length(as)))

N <- 100000
B <- sample(1:length(as),N,replace=TRUE,prob=rr)

X <- rep(0,N)
for(i in 1:N){
	X[i] <- rbeta(1,as[B[i]],bs[B[i]])
}
# 乱数から、平均と分散の推定をする
mean(X)
var(X)
# 作った関数を使ってみる
EE <- VV <- 0
for(i in 1:length(as)){
	EE <- EE + rr[i] * (as[i])/(as[i]+bs[i])
}
for(i in 1:length(as)){
	VV <- VV + rr[i] * my.beta.2moment(as[i],bs[i],EE)
}
EE
VV