如何拟合Dirichlet分布的有限混合

Question

我有一个组合样本,我想拟合Dirichlet分布的有限混合.更准确地说,请考虑以下示例:

library(gtools)
set.seed(1)
PROB = c(0.25, 0.15, 0.60)
ALPHA = list(
  c(1,1,1),
  c(2,1,1),
  c(1,1,20)
)
size = 500

N = sapply(1:3, function(i, z) sum(z == i),
           sample(1:3, size, prob = PROB, replace = TRUE))

X = do.call('rbind', 
            sapply(1:3, function(i, N) 
              rdirichlet(N[i], ALPHA[[i]]), N))[sample(1:size),]

X包含由三部分单形中定义的Dirichlet分布混合生成的样本.该混合物的第一个Dirichlet分量具有参数(1,1,1),第二个分量具有参数(2,1,1)和第三个(1,1,20).混合概率为0.25,0.15,0.60.我想从样本中检索这些参数.

你怎么会找到这个参数？

Answer 1

根据theta1 = log(p1/p3),theta2 = log(p2/p3)和所有9个alpha参数的日志进行重新参数化,然后使用带有method ="BFGS"的optim()最大化对数似然,如果使用初始值足够接近用于模拟数据的参数值.至少,Hessian的所有特征值都是负的,初始值的微小变化导致相同的最优值.

repar <- function(theta) {
  p <- exp(theta[1])
  p[2] <- exp(theta[2])
  p[3] <- 1
  p <- p/sum(p)
  alpha <- matrix(exp(theta[3:11]),3,3,byrow=TRUE)
  list(p=p,alpha=alpha)
}
logL <- function(theta,x) {
  par <- repar(theta)
  p <- par$p
  alpha <- par$alpha
  terms <- 0
  for (i in 1:length(p)) {
    terms <- terms + p[i]*ddirichlet(x,alpha[i,])
  }
  -sum(log(terms))
}
start <- c(log(c(.25,.15)/.6), log(c(1,1,1, 2,1,1, 1,1,20)))
fit <- optim(start,logL,x=X,hessian=TRUE,method="BFGS")
repar(fit$par)
eigen(fit$hessian)$val
fit2 <- optim(start+rnorm(11,sd=.2),logL,x=X,hessian=TRUE,method="BFGS")
repar(fit2$par)