I am supposed to write a function for Stirling Numbers of the Second Kind, given by the formula:
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?
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.
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")
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With