The goal
在固定的Mu上优化跨数据组的公共分散参数，其中Mu按组变化.
The problem
我有n=10组数据，我假设每组数据都是负二项分布的随机样本，固定的$Mu=mi$.
# Simulate data in hand:
# N=10 vectors of variable length
# fixed data:
N <- 10
ns <- sample(1:100, N)
# fixed means:
mus <- sample(1:100, N)
list_v <- lapply(1:N, function(i){ rnbinom(ns[i], size=3, mu=mus[i]) })
list_mu <- as.list(mus)
我想估计一个所有群体都共有的离散度(size)参数.所以我想要优化联合最大似然法，超过size个参数.我编写了一个函数negjloglik_nbinom，它可以处理变化的mu个参数:
# Define functions:
loglik_nbinom <- function(v, size, m){
log(dnbinom(v, mu=m, size=size))
}
neg_jloglik_nbinom <- function(disp, v_list, mu_list){
# Likelihood for mean m
ind_lls <- list()
for(i in 1:length(v_list)){
ind_lls[[i]] <- loglik_nbinom(size=disp, v=v_list[[i]], m=mu_list[[i]])
}
# log(product(likelihoods)) == sum(log(likelihoods))
(-1)*sum(unlist(ind_lls))
}
然后我try 将其传递给bbmle::mle2:
fit <- bbmle::mle2(minuslogl=neg_jloglik_nbinom,
start=c(disp=3),
fixed=list(v_list=list_v,
mu_list=list_mu)
)
这引发了一个奇怪的错误:
Error in bbmle::mle2(minuslogl = neg_jloglik_nbinom, start = c(disp = 3), :
some named arguments in 'fixed' are not arguments to the specified log-likelihood function:v_list1, v_list2, v_list3, v_list4, v_list5, v_list6, v_list7, v_list8, v_list9, v_list10, v_list11, v_list12, v_list13, v_list14, v_list15, v_list16, v_list17, v_list18, v_list19, v_list20, v_list21, v_list22, v_list23, v_list24, v_list25, v_list26, v_list27, v_list28, v_list29, v_list30, v_list31, v_list32, v_list33, v_list34, v_list35, v_list36, v_list37, v_list38, v_list39, v_list40, v_list41, v_list42, v_list43, v_list44, v_list45, v_list46, v_list47, v_list48, v_list49, v_list50, v_list51, v_list52, v_list53, v_list54, v_list55, v_list56, v_list57, v_list58, v_list59, v_list60, v_list61, v_list62, v_list63, v_list64, v_list65, v_list66, v_list67, v_list68, v_list69, v_list70, v_list71, v_list72, v_list73, v_list74, v_list75, v_list76, v_list77, v_list78, v_list79, v_list80, v_list81, v_list82, v_list83, v_list84, v_list85, v_list86, v_list87, v_list88, v_list89, v_list90, v_list91,
What works:
如果v_list和mu_list不是作为函数参数传递，而是neg_jloglik_nbinom在环境中找到它们，则优化最终不会成为问题.这看起来并不理想，但如果有必要，我会接受的！
# Rewrite objective function without list args:
neg_jloglik_nbinom <- function(disp){
# Likelihood for mean m
ind_lls <- list()
for(i in 1:length(v_list)){
ind_lls[[i]] <- loglik_nbinom(size=disp, v=v_list[[i]], m=mu_list[[i]])
}
# log(product(likelihoods)) == sum(log(likelihoods))
(-1)*sum(unlist(ind_lls))
}
# Assign lists to vars in environment:
v_list=list_v
mu_list=list_mu
# Compute optimization without specifying any fixed parameters:
fit <- bbmle::mle2(minuslogl=neg_jloglik_nbinom,
start=c(disp=3))