Quesiton 9, pg. 89

Stirling’s approximation provides a powerful tool for approximating the factorial of large numbers with a high degree of accuracy.

A more refined inequality for approximating \(n!\) is given by

\[ \sqrt{2\pi n}\left(\dfrac{n}{e}\right)^n e^{1/(12n+1)} < n! < \sqrt{2\pi n}\left(\dfrac{n}{e}\right)^n e^{1/(12n)}\]

Write a computer program to illustrate this inequality for n = 1 to 9

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
x <- c(1:9)
x_factorial <- factorial(x)
upper <- c()
lower <- c()

lower_fn <- function(n){
  v1 <- sqrt(2 * n * pi)
  v2 <- (n/exp(1)) ** n
  v3 <- exp(1) ** (1/(12 * n + 1))
  
  round(v1*v2*v3,4)
}

upper_fn <- function(n){
  v1 <- sqrt(2 * n * pi)
  v2 <- (n/exp(1)) ** n
  v3 <- exp(1) ** (1/(12 * n))
  
  round(v1*v2*v3, 4)
}

for (i in x){
  upper <- append(upper, upper_fn(i))
  lower <- append(lower, lower_fn(i))
}

sterlings_approx <- tibble(x, lower, x_factorial, upper)

sterlings_approx %>% mutate(ratio = lower/upper, difference = upper - lower)
## # A tibble: 9 × 6
##       x       lower x_factorial       upper        ratio     difference
##   <int>       <dbl>       <dbl>       <dbl>        <dbl>          <dbl>
## 1     1      0.9959           1      1.0023 0.9936146862 6.400000000e-3
## 2     2      1.9973           2      2.0007 0.9983005948 3.400000000e-3
## 3     3      5.9961           6      6.0006 0.9992500750 4.500000000e-3
## 4     4     23.9908          24     24.001  0.9995750177 1.020000000e-2
## 5     5    119.9699         120    120.0026 0.9997275059 3.270000000e-2
## 6     6    719.8722         720    720.0092 0.9998097247 1.370000000e-1
## 7     7   5039.3347        5040   5040.0406 0.9998599416 7.059000000e-1
## 8     8  40315.8881       40320  40320.2178 0.9998926171 4.329700000e+0
## 9     9 362850.5534      362880 362881.3779 0.9999150563 3.082450000e+1

Here I included the ratio and difference between the lower and upper bounds. We can see that the ration gets closer to 1 but the difference is getting larger between the lower and upper bounds. Running it with larger values of n we can see that the ratio continuing to get closer to 1. The accuracy of Sterling’s approximation improves as n gets larger.

library(tidyverse)

x <- seq(50,100,10)
x_factorial <- factorial(x)
upper <- c()
lower <- c()

lower_fn <- function(n){
  v1 <- sqrt(2 * n * pi)
  v2 <- (n/exp(1)) ** n
  v3 <- exp(1) ** (1/(12 * n + 1))
  
  round(v1*v2*v3,4)
}

upper_fn <- function(n){
  v1 <- sqrt(2 * n * pi)
  v2 <- (n/exp(1)) ** n
  v3 <- exp(1) ** (1/(12 * n))
  
  round(v1*v2*v3, 4)
}

for (i in x){
  upper <- append(upper, upper_fn(i))
  lower <- append(lower, lower_fn(i))
}

sterlings_approx <- tibble(x, lower, x_factorial, upper)

sterlings_approx %>% mutate(ratio = lower/upper, difference = upper - lower)
## # A tibble: 6 × 6
##       x           lower     x_factorial           upper        ratio
##   <dbl>           <dbl>           <dbl>           <dbl>        <dbl>
## 1    50 3.041400953e 64 3.041409320e 64 3.041409388e 64 0.9999972268
## 2    60 8.320971191e 81 8.320987113e 81 8.320987220e 81 0.9999980737
## 3    70 1.197855481e100 1.197857167e100 1.197857177e100 0.9999985845
## 4    80 7.156937986e118 7.156945705e118 7.156945743e118 0.9999989161
## 5    90 1.485714698e138 1.485715964e138 1.485715970e138 0.9999991435
## 6   100 9.332615095e157 9.332621544e157 9.332621570e157 0.9999993061
## # ℹ 1 more variable: difference <dbl>