Search code examples
rmath

R: Plotting a Combinatorial Function

I am working in the R programming language.

I am trying to make a plot from n=0 to n = 1000 for the following function (link: https://www.youtube.com/watch?v=iH2kATv49rc&t=7s):

enter image description here

First, I tried to define the P_00(2n) part:

combinatorial_square <- function(n) {
  (choose(2*n, n)^2) * ((1/4)^(2*n))
}

Then, I tried to write another function that performs the cumulative sum of the above function:

cumulative_sum <- function(N) {
  s <- 0
  for (n in 1:N) {
    s <- s + (choose(2*n, n)^2) * ((1/4)^(2*n))
  }
  return(s)
}

Then, I tried to plot it over a range of N's:

N <- 1:250
y <- sapply(N, cumulative_sum)

plot(N, y, type = "l", main = "Plot of cumulative_sum function", xlab = "N", ylab = "cumulative_sum(N)")

enter image description here

My Problem: I can not seem to plot this function for larger values of N (I think this might be because the computer is not able to calculate large combinatorial terms?):

 N <- 1:1000
    y <- sapply(N, cumulative_sum)

> tail(y)
[1] NaN NaN NaN NaN NaN NaN

Is there something I can do in R to approximate these larger factorials?

Currently, I was thinking of using some mathematical method to approximate the larger factorial terms involved in combinatorial expressions (e.g. https://en.wikipedia.org/wiki/Stirling%27s_approximation) - that is, re-write the cumulative_sum function so that it uses Stirling Approximation for each factorial term.

But is there an easier way to do this?

Thanks!

Solution

  • You got NAN since choose(2*n, n)^2 explodes for large n, but you should have way to circumvent that issue for big numbers.

    You can rearrange the terms in the math expression and make sure that each item is less than 1 and decay in the power law (in your case, it is power of 2), and then you will make it "converge" (not the same concept of "convergence" in series sum of the math domain), for example

    f <- Vectorize(function(n) {
  v <- seq_len(n)
  prod((n + v) / (4 * v))^2
})

    and then run (up to 10000 for exmaple)

    N <- 1:10000
P <- cumsum(f(N))
plot(N, P,  type = "l", main = "Plot of cumulative_sum function", xlab = "N", ylab = "cumulative_sum(N)")

    and you will obtain

    enter image description here