ディリクレ分布とカイ自乗分布とその歪み

• 総数を固定しなくても…

• メモ
  • は
  • (スターリングの公式)より
  • であって
  • に注意したり、,に注意したりして変形すると
  • 
• この式の移行関係には、対数尤度比がカイ自乗分布に従うという話がある

k <- 5

library(MCMCpack)
f <- c(rdirichlet(1,runif(k)))

N <- 10000
n.iter <- 1000
chi2s <- LLs <- rep(0,n.iter)
for(i in 1:n.iter){
	x <- sample(1:k,N,replace=TRUE,prob=f)
	tab <- tabulate(x,nbins = k)
	tab.exp <- N*f
	chi2 <- sum((tab-N*f)^2/(N*f))

	LL <- -2*(lgamma(N+1)-sum(lgamma(tab+1))+sum(tab * log(f)) - (lgamma(N+1)-sum(lgamma(tab.exp+1))+sum(tab.exp * log(f))))

	chi2s[i] <- chi2
	LLs[i] <- LL

}

plot(chi2s,LLs)

library(rgl)
library(MCMCpack)
R <- 3

N <- 100

f <- runif(R)
f <- f/sum(f)
#f <- rep(1/R,R)
n.iter <- 10000
X <- rdirichlet(n.iter,rep(1,R)) * N 
#X <- X + matrix(rnorm(n.iter*R),ncol=R)
v <- lgamma(apply(X,1,sum)+1) - apply(lgamma(X+1),1,sum) + apply(t(X) * log(f),2,sum)
v.2 <- apply((t(X)-f*N)^2/(f*N),2,sum)
v.2 <- v.2 - min(v.2)+0.01

d <- apply((t(X)-f*N)^2,2,sum)

plot(v,v.2)
cols <- (v * 0.5) %%3
#cols <- (d * 1.3) %%3

plot3d(X,col=cols)
open3d()
plot3d(X,col=(d * 0.01) %%3)

plot3d(X,col=(v.2/(d)) %% 3)

plot3d(X,col=gray((max(d/v.2)-d/v.2)/(max(d/v.2)-min(d/v.2))))