r - 如何拟合狄利克雷分布的有限混合

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 包含从 3 部分单纯形中定义的混合 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)