R：绘制组合函数

Question

我正在使用 R 编程语言工作。

我正在尝试为以下函数绘制从 n=0 到 n = 1000 的图（链接：https://www.youtube.com/watch?v =iH2kATv49rc&t=7s ）：

首先，我尝试定义 P_00(2n) 部分：

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

然后，我尝试编写另一个函数来执行上述函数的累积和：

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

然后，我尝试将其绘制在 N 的范围内：

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)")

我的问题：我似乎无法针对较大的 N 值绘制此函数（我认为这可能是因为计算机无法计算大的组合项？）：

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

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

我可以在 R 中做些什么来近似这些较大的阶乘吗？

目前，我正在考虑使用某种数学方法来近似组合表达式中涉及的较大阶乘项（例如https://en.wikipedia.org/wiki/Stirling%27s_approximation） - 也就是说，重写cumulative_sum函数，以便它对每个阶乘项使用斯特林近似。

但有没有更简单的方法来做到这一点呢？

谢谢！

Answer 1

诀窍是使用对数。Base R 具有beta、gamma和choose函数的对数版本factorial。这些版本以 为前缀l，在数值上是稳定的，当数字变得非常大时，它们的价值是无价的。相关函数是lchoose.

首先测试函数的结果。

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

lcombinatorial_square <- function(n) {
  2*lchoose(2*n, n) + 2*n * log(1/4)
}

N <- 1:250
y1 <- sapply(N, combinatorial_square)
y2 <- sapply(N, lcombinatorial_square)

# all values are equal to floating-point accuracy
mapply(all.equal, log(y1), y2) |> all()
#> [1] TRUE
mapply(all.equal, y1, exp(y2)) |> all()
#> [1] TRUE

# but only 1% are exactly equal
# (note that the vectors log(y1) and y2 have the same attributes
# therefore the two instructions that follow are equivalent)
mapply(identical, y1, exp(y2)) |> mean()
#> [1] 0.012
mapply(`==`, y1, exp(y2)) |> mean()
#> [1] 0.012

# (note also that if the test is on log(y1) 3% are exactly equal)
mapply(identical, log(y1), y2) |> mean()
#> [1] 0.028
mapply(`==`, log(y1), y2) |> mean()
#> [1] 0.028

# now test with larger numbers
N <- 1:1000
y2 <- sapply(N, lcombinatorial_square)
y2 |> exp() |> is.finite() |> all()
#> [1] TRUE
y2 |> exp() |> max()
#> [1] 0.25
y2 |> exp() |> sum()
#> [1] 2.265321

^{创建于 2023 年 11 月 20 日，使用reprex v2.0.2}

最后，计算总和N并绘制总和。

lcombinatorial_square <- function(n) 2*lchoose(2*n, n) + 2*n * log(1/4)

# cumulative_sum_2a <- function(N) {
#   s <- 0
#   for (n in seq.int(N)) {
#     e <- exp(lcombinatorial_square(n))
#     s <- s + e
#   }
#   s
# }

cumulative_sum_2 <- function(N) {
  N |> 
    seq.int() |>
    sapply(lcombinatorial_square) |> 
    exp() |> 
    sum()
}

N <- 1:1000
y <- sapply(N, cumulative_sum_2)
plot(N, y, type = "l", main = "Plot of cumulative_sum function", 
     xlab = "N", ylab = "cumulative_sum(N)")

^{创建于 2023 年 11 月 20 日，使用reprex v2.0.2}