使用具有多个参数和索引的lapply / sapply减少函数内的嵌套循环

Posted 2021-03-27

我有两个观察向量，我想要计算哪些参数最适合这些观察。因此，我已经定义了似然函数。我正在努力从嵌套for循环到使用sapply。我的目标是返回观察“芽”概率的（50x1）向量。

> bud
 [1] 1.72 1.12 1.01 1.14 1.56 1.04 1.00 1.39 1.38
[10] 0.99 1.66 1.55 2.26 0.86 1.45 0.90 1.58 1.09
[19] 0.97 0.98 1.16 1.20 1.17 1.59 1.00 0.92 1.24
[28] 0.98 1.29 1.36 1.80 1.42 1.67 1.13 0.76 0.93
[37] 1.03 1.29 1.09 1.64 1.27 1.14 1.60 1.00 0.91
[46] 1.40 1.15 2.03 1.34 1.37
> deltakere
 [1] 8 5 7 6 9 5 5 9 9 6 5 8 8 6 6 5 7 6 7 7 9 6 6 8
[25] 8 6 6 5 9 9 6 7 9 5 5 6 6 8 7 9 8 6 9 8 5 8 8 6
[49] 9 7

当前代码（尽管工作缓慢）：

# first equation in my likelihood function
eqOne <- function(i,theta1, theta2, lambda){
        vektorN = (1:100)
        loopVerdi = 0
            if(deltakere[i]>=2){        
                for (N in deltakere[i]:20000) {
                    density = ((N)*(N-1)*(pWEI2(bud[i], theta1, theta2))^(N-2) *
                        pWEI2(bud[i], theta1, theta2, lower.tail = FALSE) *
                        dWEI2(bud[i], theta1, theta2)) /
                        (1-(pWEI2(r, theta1, theta2))^N)
                    ngittN =
                        dbinom(deltakere[i], size = N, 
                               prob = pWEI2(r, theta1, theta2, lower.tail = FALSE))
                    sN = 
                        dpois(N, lambda)
                    loopVerdi = loopVerdi + (density * ngittN * sN)}
            return(loopVerdi)}}

# here is the second equation in the likelihood estimate.
eqTwo <- function(theta1, theta2, lambda) {
        firstPart <- 1:length(bud)
            for (i in 1:length(bud)){
                firstPart[i] <- eqOne(i,theta1, theta2, lambda)}
        return(firstPart)}

Vector <- eqTwo(1,2,7)
> Vector
 [1] 0.071126958 0.174198273 0.125381705 0.202950594
 [5] 0.116409283 0.156245797 0.143071864 0.156875579
 [9] 0.157670857 0.139149430 0.052085368 0.138248017
[13] 0.001884680 0.073326986 0.161989881 0.103352207
[17] 0.125124465 0.186201434 0.103156346 0.108636007
[21] 0.116075468 0.214339615 0.209921459 0.121433257
[25] 0.084230719 0.102485943 0.216127393 0.135710730
[29] 0.153834163 0.158552633 0.036673600 0.194744214
[33] 0.079498794 0.175511921 0.049722001 0.107646750
[37] 0.159465393 0.198467993 0.169062724 0.089375280
[41] 0.196173239 0.202950594 0.102923165 0.084230719
[45] 0.107521392 0.190136547 0.159243477 0.008274964
[49] 0.158457042 0.210361176

这段代码产生一个（50x1）向量，具有观察“芽”值的概率，这就是我想要的。

这是我尝试从for循环到sapply方法

funDeltakere <- function(i, theta1, theta2, lambda){ 
    function(N, theta1, theta2, lambda){
                    density = ((N) *(N-1)*(pWEI2(bud[i], theta1, theta2))^(N-2) *
                        pWEI2(bud[i], theta1, theta2, lower.tail = FALSE) *
                        dWEI2(bud[i], theta1, theta2)) /
                        (1-(pWEI2(r, theta1, theta2))^N)
                    ngittN =
                        dbinom(deltakere[i], size = N, 
                               prob = pWEI2(r, theta1, theta2, lower.tail = FALSE))
                    sN = dpois(N, lambda)
                    return(density * ngittN * sN)}}

eqOne <- function(i,theta1, theta2, lambda){
            if(deltakere[i]>=2){        
                loopVerdi <- sapply(deltakere[i]:20000, funDeltakere(i, theta1, theta2, lambda))
            return(loopVerdi)}} 

eqTwo <- function(theta1, theta2, lambda) {
            firstPart <- 1:length(bud)
                for (i in 1:length(bud)){
                    firstPart[i] <- eqOne(i,theta1, theta2, lambda)
                }
            return(firstPart)
            }
Vector <- eqTwo(1,2,7)

pWEI2（bud [i]，theta1，theta2）出错：参数“theta1”缺失，没有默认值

看起来这些参数并没有在我的代码中通过我的第一个函数。任何有关如何确保我的参数通过所有功能的提示将不胜感激！

答案

我设法解决了我的问题。

bud <- runif(10, 1.1,1.55)
deltakere <- round(runif(10,5,9))
r=0.5

foo_inner <- function(i, theta1, theta2, lambda){
    function(N){
    density = (
        (N) *
        (N-1) * 
        (pWEI2(bud[i], theta1, theta2))^(N-2) *
        pWEI2(bud[i], theta1, theta2, lower.tail = FALSE) *
        dWEI2(bud[i], theta1, theta2)) /
        (1-(pWEI2(r, theta1, theta2))^N)

    ngittN =
        dbinom(deltakere[i], size = N, 
               prob = pWEI2(r, theta1, theta2, lower.tail = FALSE))
    sN = 
        dpois(N, lambda)

    return(density * ngittN * sN)
    }}

foo_outer <- function(theta1, theta2, lambda){
            function(i){
            listeObs <- sapply(deltakere[i]:20000, foo_inner(i, theta1, theta2, lambda))
            return(sum(listeObs))
            }}

eqThree <- function(theta1, theta2, lambda){
    secondPart <- pbsapply(1:length(bud), foo_outer(theta1, theta2, lambda))

    LL <- -sum(log(secondPart))
    return(LL)
}

result_mle <- mle(minuslogl = eqThree, start=list(theta1 = 1,
                                                  theta2 = 2,
                                                  lambda = 7),
                  method="L-BFGS-B", lower=c(0.1,0.1,5),
                  nobs = length(bud))