2017-07-14

ディリクレ過程で二項分布その２

昨日の記事でRのDPpackageのDPbetabinom()関数の使い方をさらった
観測データが320のオブジェクトについて9回繰り返し試行をしているわけだが、9回の繰り返し試行で、それの「表確率」の推定はそれほど精度がよくできるものではない。
個々のオブジェクトでの試行回数が大きくなれば、そのオブジェクトの成功確率の推定値精度は上がる。
その条件を満足したまま、たくさんのオブジェクトで観察すれば、個々のオブジェクトの成功確率を推定した上で、その点推定値の分布を $G_0$ とみなして差支えないだろう。
逆に、オブジェクト数が少なければ、すかすかのデータなのでどのような $G_0$ かの推定精度は下がらざるを得ない。
一方、各オブジェクトの試行回数が少ないときは、個々のオブジェクトの成功確率の推定幅は広い。それは、複数のオブジェクトが同一クラスタに帰属すると考えてよい可能性を上げる。
オブジェクト数が多いときには、クラスタ数が少なくてもそれなりに説明できることになる。オブジェクト数が少ないときにもクラスタ数が少なくてもそれなりに説明できる。両者の違いは、オブジェクト数が多い方が、クラスタごとの成功確率推定値の精度がよくなること(たぶん)と、いくつのクラスタにするべきかの精度が上がること(たぶん)
Rでやってみる。オブジェクト数を10,100,1000。各オブジェクトの試行回数を10,100,10000として、 $\alpha$ の事前分布をガンマ分布(a0=1,b0=1)とした場合を示す。
「真実」はクラスタ数５で、その割合と、各クラスタの成功確率を示した上で、DPbetabinom()の結果のp分布を示す。

library(DPpackage)
library(MCMCpack)

################ 背景分布
n.cl <- 5 # クラスタ数
rs <- rdirichlet(1,rep(5,n.cl)) # クラスタごとの割合
ps <- sort(runif(n.cl)*0.4+0.2) # 各クラスタの表確率

plot(ps,rs,type="h",xlim=c(0,1))
################
################推定条件
prior4<-list(a0=1,
						b0=1,
						a1=1,
						b1=1)
# Initial state
state <- NULL
# MCMC parameters
ngrid <- 100
nsave <- 1000
mcmc2 <- list(nburn=500,
						 nsave=nsave,
						 nskip=3,
						 ndisplay=100,
						 ngrid=ngrid)
################
################観察条件

N.samples <- c(10,100,1000)
N.trials <- c(10,100,10000)

NN.st <- expand.grid(N.samples,N.trials)
###############

###############推定結果記録
densp.ms <- matrix(0,length(NN.st[,1]),ngrid)
cl.states <- matrix(0,length(NN.st[,1]),nsave)
alphas <- matrix(0,length(NN.st[,1]),nsave)
###############
for(ii in 1:length(NN.st[,1])){
	N.sample <- NN.st[ii,1] # サンプル数
	N.trial <- NN.st[ii,2] # 各サンプルでの試行回数

	# クラスタを割り当てる
	cl <- sample(1:n.cl,N.sample,replace=TRUE,prob=rs)
	# 観察データを作る
	y <- matrix(0,N.sample,2)
	for(i in 1:N.sample){
		tmp <- sample(1:0,N.trial,replace=TRUE,prob=c(ps[cl[i]],1-ps[cl[i]]))
		y[i,1] <- sum(tmp)
		y[i,2] <- N.trial
	}


	# Fitting the model

	fit4 <- DPbetabinom(y=y,ngrid=100,
	                       prior=prior4,
	                       mcmc=mcmc2,
	                       state=state,
	                       status=TRUE)

	densp.ms[ii,] <- fit4$densp.m
	cl.states[ii,] <- fit4$save.state$thetasave[,1]
	alphas[ii,] <- fit4$save.state$thetasave[,2]

}

par(mfrow=c(3,3))
for(ii in 1:length(NN.st[,1])){
	ttl <- paste("Nsample",NN.st[ii,1]," Niter ",NN.st[ii,2])
	plot(fit4$grid,densp.ms[ii,],type="l",main=ttl)
}

条件を変えて、 $\alpha=100$ のようにクラスタ数がおおくなるように指定すると以下のようになる。実際、推定クラス多数は、真実の5よりも大幅に大きくななる

prior3<-list(alpha=100,
						a1=1,
						b1=1)
for(ii in 1:length(NN.st[,1])){
	ttl <- paste("Nsample",NN.st[ii,1]," Niter ",NN.st[ii,2])
	hist(cl.states[ii,],main=ttl)
}