Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to plot user-defined functions in R?

I am supposed to write a function for Stirling Numbers of the Second Kind, given by the formula:

enter image description here

For this, I have written the following function in R:

stirling <- function(n, k)
{
  sum = 0
  for (i in 0:k)
  {
    sum = sum + (-1)^(k - i) * choose(k, i) * i^n
  }
  sum = sum / factorial(k)
  return(sum)
}

The next part of the question is to "create a plot for n = 20, k = 1,2,...,10". I did some research and I think the methods curve or plot might help me. However, I am guessing these methods are used when y is of the form f(x) (i.e. a single argument). But here, I have two arguments (n and k) in my function stirling so I am not sure how to approach this.

Also, I tried converting the values of k (0, 1, 2..., 10) to a vector and then passing them to stirling, but stirling won't accept vectors as input. I am not sure how to modify the code to make stirling accept vectors.

Any suggestions?

like image 610
lebowski Avatar asked Sep 03 '16 03:09

lebowski


People also ask

How do I plot multiple functions in R?

To draw multiple curves in one plot, different functions are created separately and the curve() function is called repeatedly for each curve function. The call for every other curve() function except for the first one should have added an attribute set to TRUE so that multiple curves can be added to the same plot.


1 Answers

Vectorize

As pointed out in the comments, you can vectorize to do this:

Vectorize creates a function wrapper that vectorizes the action of its argument FUN. Vectorize(FUN, vectorize.args = arg.names, SIMPLIFY = TRUE, USE.NAMES = TRUE)

(vstirling <- Vectorize(stirling))
# function (n, k) 
# {
# args <- lapply(as.list(match.call())[-1L], eval, parent.frame())
# names <- if (is.null(names(args))) 
#     character(length(args))
# else names(args)
# dovec <- names %in% vectorize.args
# do.call("mapply", c(FUN = FUN, args[dovec], MoreArgs = list(args[!dovec]), 
#    SIMPLIFY = SIMPLIFY, USE.NAMES = USE.NAMES))
# }

so vstirling() is the vectorized version of stirling().

vstirling(20, 1:10)
 # [1] 1.000000e+00 5.242870e+05 5.806064e+08 4.523212e+10 7.492061e+11 4.306079e+12 1.114355e+13 1.517093e+13
 # [9] 1.201128e+13 5.917585e+12

Now all that is left is creating a plot:

plot(x = 1:10, y = vstirling(20, 1:10), ylab = "S(20, x)", xlab = "x")
like image 117
KenHBS Avatar answered Sep 23 '22 08:09

KenHBS