Добавление функции в цикл for в R

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), чей первый элемент является mean используется в dnorm части, а остальные используются в 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 

Согласно комментарий Николы, да, это медленнее.

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

Около 4-5 раз медленнее, т.е. ~ 5мс против ~ 1мс работать с 1000 средствами.

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 

Edit: Были обновлены до vectorise код, добавлен тест