如何拟合Dirichlet分布的有限混合
如何拟合Dirichlet分布的有限混合
提问于 2016-05-23 08:50:07
我有一个组合样本，我想拟合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。我想从示例中检索这些参数。
你怎么找到这个参数？
回答 1
Stack Overflow用户
回答已采纳
发布于 2016-06-05 19:58:19
根据theta1=log(p1/p3)、theta2=log(p2/p3)和所有9个alpha参数的日志重新参数化，然后使用optim()和method="BFGS“最大限度地提高日志的可能性，如果使用与模拟数据所使用的参数值足够接近的初始值，则似乎是可行的。至少，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)页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/37386550
复制相关文章
相似问题
腾讯云开发者
