From bcc9220ee32382601561e25de7610fab20752745 Mon Sep 17 00:00:00 2001
From: "Hyeong Kyun (Daniel) Park" <park1119@purdue.edu>
Date: Tue, 2 Nov 2021 23:56:06 -0400
Subject: [PATCH] finished chapter 6

---
 ch6.R | 169 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 169 insertions(+)
 create mode 100644 ch6.R

diff --git a/ch6.R b/ch6.R
new file mode 100644
index 0000000..d7a13ec
--- /dev/null
+++ b/ch6.R
@@ -0,0 +1,169 @@
+#############
+# Chapter 6.2 (Probability Inequality)
+
+## Compare Chebyshev's and Chernoff's bounds
+
+# R code to compare the probability bounds
+
+library(pracma)
+
+epsilon <- 0.1
+sigma <- 1;
+N <- logspace(1,3.9,50)
+p_exact <- 1-pnorm(N**(1/2)*epsilon/sigma, 0, 1)
+p_cheby <- sigma**2. / (epsilon**2*N)
+p_chern <- exp(-epsilon**2*N/(2*sigma**2))
+
+plot(log(N), log(p_exact), pch=1, col="orange", lwd=2, xlab="log(N)", ylab="log(Probability)")
+points(log(N), log(p_cheby), pch=15, col="green", lwd=2)
+lines(log(N), log(p_chern), pch=19, col="blue", lwd=2)
+legend("bottomleft", c("Exact","Chebyshev","Chernoff"), fill=c("orange", "green", "blue"))
+
+#############
+# Chapter 6.3 (Law of Large Numbers)
+
+## Weak law of large numbers
+
+# R code to illustrate the weak law of large numbers
+
+library(pracma)
+
+p <- 0.5
+Nset <- as.integer(round(logspace(2,5,100)))
+x <- matrix(rep(0, 1000*length(Nset)), nrow=1000)
+for (i in 1:length(Nset)) {
+  N = Nset[i]
+  x[,i] <- rbinom(1000, N, p) / N
+}
+Nset_grid <- repmat(Nset, m=1, n=1000)
+
+semilogx(Nset_grid, x, col='black', pch=19)
+points(Nset, p + 3*(((p*(1-p))/Nset)^(1/2)), col='red', pch=19, lwd=1)
+points(Nset, p - 3*(((p*(1-p))/Nset)^(1/2)), col='red', pch=19, lwd=1)
+
+#############
+# Chapter 6.4 (Central Limit Theorem)
+
+## PDF of the sum of two Gaussians
+
+# Plot the PDF of the sum of two Gaussians
+
+library(pracma)
+
+n <- 10000
+K <- 2
+Z <- rep(0, n)
+
+for (i in 1:K) {
+  X <- runif(n, min=1, max=6)
+  Z <- Z + X
+}
+hist(Z,breaks=(K-0.5):(6*K+0.5),freq=FALSE)
+
+# Visualize convergence in distribution
+
+library(pracma)
+
+N <- 10
+N <- 50
+x <- linspace(0, N, 1001)
+p <- 0.5
+p_b <- dbinom(x, N, p)
+p_n <- dnorm(x, N*p, (N*p*(1-p))**(1/2))
+
+c_b <- pbinom(x, N, p)
+c_n <- pnorm(x, N*p, (N*p*(1-p))**(1/2))
+
+plot(x, p_n, lwd=1, col='red')
+lines(x, p_b, lwd=2, col='black')
+legend("topright", c('Binomial', 'Gaussian'), fill=c('black', 'red'))
+
+# Poisson to Gaussian: convergence in distribution
+
+library(pracma)
+
+N <- 4
+# N = 10
+# N = 50
+
+x <- linspace(0,2*N,1001)
+lambda <- 1
+p_b <- dpois(x, N*lambda)
+p_n <- dnorm(x, N*lambda, sqrt(N*lambda))
+
+c_b <- ppois(x, N*lambda);
+c_n = pnorm(x, N*lambda, sqrt(N*lambda));
+
+plot(x, p_n, col="red")
+lines(x, p_b, col="black")
+legend("topright", c('Poisson', 'Gaussian'), fill=c('black', 'red'))
+
+plot(x, c_n, col="red")
+lines(x, c_b, col="black")
+legend("topright", c('Poisson', 'Gaussian'), fill=c('black', 'red'))
+
+# Visualize the Central Limit Theorem
+
+library(pracma)
+
+N <- 10
+x <- linspace(0,N,1001)
+p <- 0.5
+p_b <- dbinom(x, N, p);
+p_n <- dnorm(x, N*p, sqrt(N*p*(1-p)));
+
+c_b <- pbinom(x, N, p);
+c_n <- pnorm(x, N*p, sqrt(N*p*(1-p)));
+
+x2 <- linspace(5-2.5,5+2.5,1001);
+q2 <- dnorm(x2,N*p, sqrt(N*p*(1-p)));
+
+plot(x, p_n, col="red")
+points(x, p_b, col="black", pch=19)
+
+polygon(c(min(x2), x2, max(x2)), c(0, q2, 0), col='lightblue')
+
+# How moment generating of Gaussian approximates in CLT
+
+library(pracma)
+
+p <- 0.5
+s <- linspace(-10,10,1001)
+MX <- 1-p+p*exp(s)
+N <- 2
+semilogy(s, (1-p+p*exp(s/N))**N, lwd=4, col="lightblue", xlim=c(-10,10), ylim=c(10**-2, 10**5))
+mu <- p
+sigma <- sqrt(p*(1-p)/N);
+MZ <- exp(mu*s + sigma^2*s**2/2);
+lines(s, MZ, lwd=4);
+legend("topleft", c('Binomial MGF', 'Gaussian MGF'), fill=c('lightblue', 'black'))
+
+
+# Failure of Central Limit Theorem at tails
+
+library(pracma)
+
+x <- linspace(-1,5,1001)
+lambda <- 1
+
+N <- 1
+f1 <- (N**(1/2)/lambda)*dgamma((x+sqrt(N))/(lambda/sqrt(N)), N, lambda)
+semilogy(x, f1, lwd=1, col='lightgray', xlim=c(-1,5), ylim=c(10**-6, 1))
+
+N <- 10
+f1 <- (N**(1/2)/lambda)*dgamma((x+sqrt(N))/(lambda/sqrt(N)), N, lambda)
+lines(x, f1, lwd=2, col='gray')
+
+N <- 100
+f1 <- (N**(1/2)/lambda)*dgamma((x+sqrt(N))/(lambda/sqrt(N)), N, lambda)
+lines(x, f1, lwd=2, col='darkgray')
+
+N <- 1000
+f1 <- (N**(1/2)/lambda)*dgamma((x+sqrt(N))/(lambda/sqrt(N)), N, lambda)
+lines(x, f1, lwd=2, col='black')
+
+g <- dnorm(x,0,1)
+lines(x, g, lwd=2, pch=1, col='red')
+
+legend("bottomleft", c('N=1', 'N=10', 'N=100', 'N=1000', 'Gaussian'), fill=c('lightgray', 'gray', 'darkgray', 'black', 'red'))
+