少し改良

• カテゴリ数を増やしつつ、その「未観測カテゴリ」の確率をプールして計算するだけなら
• をカテゴリ数仮説の相対尤度とし、「未観測カテゴリ１個分」の生起確率期待値はなので、の未観測カテゴリ数分を合算すればについて、カテゴリ数仮説の相対尤度重みづけ平均を出し、その余りを観測度数+1に応じて振り分ければよい
• 仮説の相対確率・尤度の総和や重みづけの計算では、仮説間の対数尤度の開きが大きいと計算がオーバーフローするので、どのみち小さい相対尤度の仮説は0になるように相対対数尤度を最大のそれがハンドリングできるように調整するのもこつ

my.nkK <- function(n,k,K){
	g1 <- lgamma(n+K)
	g2 <- lgamma(2*n+K)
	a <- K-k
	b <- n+K
	return(list(G = g1-g2,ab = a/b))
}

my.nkK.series <- function(n,k,maxK){
	K <- (k):maxK
	Gs <- abs <- rep(0,length(K))
	for(i in 1:length(K)){
		tmp <- my.nkK(n,k,K[i])
		Gs[i] <- tmp$G
		abs[i] <- tmp$ab
	}
	Gs <- Gs-max(Gs)+log(.Machine$integer.max-10)
	#Gs <- Gs-max(Gs)+15
	#pr <- exp(Gs)
	#pr. <- pr/sum(pr)
	ret <- sum(exp(Gs+log(abs)))/sum(exp(Gs))
	#ret <- sum(pr.*abs)
	return(list(pr = ret, Gs=Gs,abs=abs))
}


my.full.dirichlet <- function(v,maxK){
	n <- sum(v)
	k <- length(v)
	out <- my.nkK.series(n,k,maxK)
	PR <- 1-out$pr
	p <- (v+1)/sum(v+1) * PR
	return(c(p,out$pr))
}

• 前の記事のソースとくらべてみれば

v <- c(7,4,3,2)
maxK <- 10000
my.full.dirichlet(v,maxK)
my.exp.dirichlet.multi(v,maxK)$expected.prob.obs.unobs

• で出力は同じで、新しい方が格段に速い