在R中的for循环中附加一个函数

Question

我有积分函数，它是几个累积概率函数和密度函数的乘积。

对于两个事件，积分函数只是累积概率和密度函数的乘积：

function(value) pnorm(value,mean = mean2 ,sd = sigma, lower.tail = TRUE)*
dnorm(value,mean = mean1, sd = sigma)

对于每个新事件，我都需要乘以另一个累积概率函数。因此，对于三个备选方案，函数变为：

function(value) pnorm(value,mean = mean2 ,sd = sigma, lower.tail = TRUE)*
 pnorm(value,mean = mean1,sd = sigma, lower.tail = TRUE)*
 dnorm(value,mean = mean0, sd = sigma)

四个人：

 function(value) pnorm(value,mean = mean3 ,sd = sigma, lower.tail = TRUE)*
 pnorm(value,mean = mean2,sd = sigma, lower.tail = TRUE)*
 pnorm(value,mean = mean1,sd = sigma, lower.tail = TRUE)*
 dnorm(value,mean = mean0, sd = sigma)

等等……

我正在尝试构建一个循环，为任意数量的事件动态创建此函数。我尝试了不同的方法来概括该功能，但到目前为止对我来说没有任何效果。关于我应该如何进行的任何想法？

Answer 1

使用我提出的方法，您将一个向量 (meanvec) 传递给函数，其第一个元素是在 dnorm 部分中使用的 mean，其他元素在 pnorm 中使用。我假设sigma 总是相同的。使用这种方法，您可以在 meanvec 参数中传递任意数量的元素。

myfun<-function(value,meanvec,sigma) {
    valuelong<-rep(value,each=length(meanvec)-1)
    ppart<-apply(matrix(pnorm(valuelong,meanvec[-1],sigma),nrow=length(meanvec)-1),2,prod)
    if (length(meanvec)>1) ppart*dnorm(value,meanvec[1],sigma) else dnorm(value,meanvec[1],sigma)
}

例子：

mean0<-1
mean1<-2
mean2<-3
mean3<-4
value<-runif(100)
sigma<-2
#here define your not generalized function with 3 pnorm
oldfun<-function(value) pnorm(value,mean = mean3 ,sd = sigma, lower.tail = TRUE)*
pnorm(value,mean = mean2,sd = sigma, lower.tail = TRUE)*
pnorm(value,mean = mean1,sd = sigma, lower.tail = TRUE)*
dnorm(value,mean = mean0, sd = sigma)
all.equal(oldfun(value),myfun(value,c(mean0,mean1,mean2,mean3),sigma))
#[1] TRUE

Answer 2

这是另一种我认为很容易遵循的方法。

这个的关键是Reduce 函数。即Reduce('+',1:3)与1+2+3相同，Reduce('*',1:3)与1*2*3相同。

do_it <- function(value, means, sigma) {

  # get the pnorm values (and then ignore the first result)
  pn_result <- pnorm(value, means, sigma, lower.tail = TRUE)[-1]

  # now get the dnorm (i.e. the first mean)
  dn_result  <- dnorm(value[1], means[1], sd = sigma[1])

  # Use reduce function to multiply all values together
  Reduce("*", c(pn_result, dn_result))

}

现在像这样使用函数：

> do_it(value = 7,   means = 2:4, sigma = 2)
[1] 0.007992577
> do_it(value = 7,   means = 2:4, sigma = 1:3)
[1] 1.222387e-06
> do_it(value = 7:9, means = 2:4, sigma = 1:3)
[1] 1.406878e-06

根据 Nicola 的评论，是的，这比较慢。

microbenchmark(
  do_it(7, 2:1000, 2),
  myfun(7, 2:1000, 2),
  time = 10000,
  unit = 'eps'
)

以 1000 次均值运行大约慢 4-5 倍，即 ~5ms 与 ~1ms。

Unit: evaluations per second
                expr         min           lq        mean       median           uq          max neval
 do_it(7, 2:1000, 2)    1006.555     2234.863     2299.33     2373.921     2409.485 2.554957e+03   100
 myfun(7, 2:1000, 2)    5627.335     9837.340    10040.78    10169.636    10497.424 1.155161e+04   100
                time 9523809.524 30776515.152 48510649.75 38461538.462 57189542.484 1.666667e+08   100

编辑：已更新以矢量化代码，添加了基准