R中生产内部求和

Summation inside summation inside production in R

我在编写要优化的函数时遇到问题,其中有两个求和和一个产生式,都具有不同的索引。为了简单起见,我将代码分成两个函数。

在第一个函数中,j 从 0 到 k:

w = function(n,k,gam){
  j = 0:k
  w = (1 / factorial(k)) * n * sum(choose(k, j * gam))
  return(w)}

在第二个函数中,k 从 0 变为 n(即固定为 10);相反,生产从 1 到 length(x):

    f = function(gam,del){
   x = mydata #vector of 500 elements
   n = 10    
 k = 0:10  
 for (i in 0:10)
        pdf = prod( sum( w(n, k[i], gam) * (1 / del + (n/x)^(n+1))
return(-pdf)}

当我尝试该函数时出现以下错误:

Error in 0:k : argument of length 0

编辑:这就是我要编写的代码

我想最大化 L(d,g) 使用 optim 和:

n固定为特定值。

解决方案

for (i in 0:10) 更改为 for ( i in 1:11 )。注意:当我复制 运行 你的代码时,我还注意到一些不相关的 bracket/parentheses 遗漏,你可能还需要修复。

说明

您的问题是 R 使用基于 1 的索引系统,而不是像许多其他编程语言或某些数学公式那样使用基于 0 的索引系统。如果你 运行 下面的代码你会得到同样的错误,它指出了问题所在:

k = 0:10
for ( i in 0:10 ) {
    print(0:k[i])
}

Error in 0:k[i] : argument of length 0

第一次迭代时出现错误,因为 k 没有 0 元素。将其与以下循环进行比较:

k = 0:10
for ( i in 1:11 ) {
    print(0:k[i])
}

[1] 0
[1] 0 1
[1] 0 1 2
[1] 0 1 2 3
[1] 0 1 2 3 4
[1] 0 1 2 3 4 5
[1] 0 1 2 3 4 5 6
[1] 0 1 2 3 4 5 6 7
[1] 0 1 2 3 4 5 6 7 8
 [1] 0 1 2 3 4 5 6 7 8 9
 [1]  0  1  2  3  4  5  6  7  8  9 10

更新

您对答案的评论阐明了您需要的一些额外信息:

Just to full understand everything, how do I know in a situation like this that R is indexing the production on x and the summation on k?

简短的回答是,这取决于您嵌套循环和函数调用的方式。更详细:

  • 当您调用 f() 时,您会在 k 的元素上启动 for 循环,因此 R 正在索引 [=19] 中的代码块=]循环(下面我重新格式化的f()版本中大括号中的所有内容)"on"k。对于 k 中的每个元素,您将 prod(...) 分配给 pdf(旁注:我不知道为什么您要在每次迭代中重写 pdf循环)
  • sum( w(n, k[i], gam) * gamma(1 / del + k[i]) * s^(n + 1)) 产生一个长度为 max(length(w(n, k[i], gam)), length(s)) 的向量(边注:注意回收!——参见 Section 2.2 of "An Introduction to R"); prod(sum( w(n, k[i], gam) * gamma(1 / del + k[i]) * s^(n + 1))) 有效索引该向量的元素
  • w(n, k[i], gam) * gamma(1 / del + k[i]) * s^(n + 1) 产生一个长度为 max(length(w(n, k[i], gam)), length(s)) 的向量; sum( w(n, k[i], gam) * gamma(1 / del + k[i]) * s^(n + 1)) 有效索引该向量的元素
  • 等等

您通过矢量化操作显式或隐式索引的内容取决于您所讨论的嵌套循环或函数调用的级别。您可能需要仔细考虑和计划何时要索引什么,这将告诉您需要如何嵌套内容。将索引变化最快的操作放在最里面的调用中。例如,在 effect 中,prod(1:3 + sum(1:3)) 将首先对 sum(1:3) 进行索引以生成该和,然后对 1:3 + sum(1:3) 进行索引以生成产品。即,sum(1:3) = 1 + 2 + 3 = 6,然后 prod(1:3 + sum(1:3)) = (1 + 6) * (2 + 6) * (3 + 6) = 7 * 8 * 9 = 504。它是就像括号在数学中的作用一样。

此外,另一方面,我不会像您在 f() 中那样从函数中引用全局变量——我在下面的代码中突出显示了您这样做的地方,并提供了一个替代方案那不行。

f = function(gam, del){
    x = mydata # don't refer to a global variable "mydata", make it an argument
    n = 10  
    s = n / x  
    k = 1:11
    for (i in 1:11){
        pdf = prod( sum( w(n, k[i], gam) * gamma(1 / del + k[i]) * s^(n + 1)))
    }
    return(-pdf)
}

# Do this instead
# (though there are still other things to fix,
# like re-writing over "pdf" eleven times and only using the last value)

f = function(gam, del, x, n = 10) {
    s = n / x
    s = n / x  
    k = 0:10
    for (i in 1:11){
        pdf = prod( sum( w(n, k[i], gam) * gamma(1 / del + k[i]) * s^(n + 1)))
    }
    return(-pdf)
}