Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Vectorized R function

I have to write a vectorized R function f that takes a vector x= (x_1, . . . , x_m) and a natural number n, and returns the value f_n(x) given by:

enter image description here

Example:

> x = seq(-1, 3, by = 0.4)
> f(x,6)   # here n=6

[1] 0.000000e+00 0.000000e+00 0.000000e+00
[4] 2.666667e-06 6.480000e-04 8.333333e-03
[7] 4.430667e-02 1.410800e-01 3.050933e-01
[10] 4.755467e-01 5.500000e-01

This is what I got:

f = function(x, n){
  s = 0
  for(j in 0:x)
    s = s + (-1)^j*choose(n, j)*(x-j)^(n-1)
  s/factorial(n-1)
}

x = seq(-1, 3, by = 0.4)
f(x,6)

Warning in 0:x: numerical expression has 11 elements: only the first used
[1] -8.333333e-03 -6.480000e-04 -2.666667e-06 2.666667e-06 6.480000e-04
[6] 8.333333e-03 4.481867e-02 1.574640e-01 4.294693e-01 9.901147e-01
[11] 2.025000e+00

Clearly it is not what it should be in the example. What did I do wrong here? TIA

EDIT: Maybe using outer and apply might help with x?


2 Answers

This is a slightly different way of doing this solely based on base R:

x = seq(-1, 3, by = 0.4)
n <- 6

fn <- function(x, n) {
  x[x <= 0] <- 0
  sapply(x, function(x) {
    Reduce(function(a, b) {
      a + (-1) ^ b * (factorial(n)/(factorial(b) * factorial(n-b))) * (x - b) ^ (n-1)
    }, seq(0, x), init = 0) * (1/factorial(n-1))
  })
}

fn(x, 6)

 [1] 0.000000e+00 0.000000e+00 0.000000e+00 2.666667e-06 6.480000e-04 8.333333e-03 4.430667e-02
 [8] 1.410800e-01 3.050933e-01 4.755467e-01 5.500000e-01
like image 52
Anoushiravan R Avatar answered Sep 19 '26 05:09

Anoushiravan R


Try this code. It can be modified to become tidier but maybe it can solve your problem in its current form. I used both base R and purr functions for iteration instead of for loop but maybe for loop alone can do the job.

library(tidyverse)
n <- 6
x <- seq(-1, 3, by = 0.4)
x[x<= 0] <- 0
seq_fun <- function(x) seq(0, x)
d <- sapply(x, seq_fun)
fun <- function(r, t) {
   sum((-1) ^ r *choose(n, r)*(r-t)^(n-1)) / factorial(n - 1)
}
as_vector(map2(d, x, fun))
like image 23
ETeddy Avatar answered Sep 19 '26 06:09

ETeddy



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!