2

我有一个组成样本，我想拟合 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 分量具有参数 (1,1,1)，第二个分量具有参数 (2,1,1)，第三个分量具有参数 (1,1,20)。混合概率为 0.25、0.15、0.60。我想从示例中检索这些参数。

你如何找到这个参数？

4

1 回答 1

6

根据 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)
于 2016-06-05T19:58:19.653 回答