f <- function(x){
  exp(-x^2)
}

integrate(f, -Inf, Inf)

sqrt(pi)

# Sums of binomials is binomial

B <- 1000

x <- rbinom(B, size = 20, prob = 0.55)
y <- rbinom(B, size = 30, prob = 0.55)
z <- x + y
par(mfrow = c(3, 1))
hist(x,  probability = TRUE, xlim = c(0, 50))
hist(y,  probability = TRUE, xlim = c(0, 50))
hist(z,  probability = TRUE, xlim = c(0, 50), main = "Sum of binomials")
# Mean and variance of the sum  
mean(z)
# Sum of the means
20 * 0.55 + 30 * 0.55
var(z)
# Sum of the variances
20 * 0.55 * 0.45 + 30 * 0.55 * 0.45

# Sums of Exponential is Gamma

x <- rexp(B, rate = 0.3)
y <- rexp(B, rate = 0.2)
z <- x + y
par(mfrow = c(3, 1))
hist(x,  probability = TRUE, xlim = c(0, 60))
hist(y,  probability = TRUE, xlim = c(0, 60))
hist(z,  probability = TRUE, xlim = c(0, 60), main = "Sum of exponentials")
# Mean and variance of the sum  
mean(z)
# Sum of the means
1/0.1 + 1/0.1
var(z)
# Sum of the variances
1/0.1^2 + 1/0.1^2
# Compare to Gamma distribution
curve(dgamma(x, shape = 2, rate = 0.1), add = TRUE, col = "red", lwd = 2)


# Sums of Normals is Normal

x <- rnorm(B, mean = 10, sd = 3)
y <- rnorm(B, mean = 20, sd = 4)
z <- x + y
par(mfrow = c(3, 1))
hist(x,  probability = TRUE, xlim = c(0, 50))
hist(y,  probability = TRUE, xlim = c(0, 50))
hist(z,  probability = TRUE, xlim = c(0, 50), main = "Sum of normals")
# Mean and variance of the sum
mean(z)
# Sum of the means
10 + 20
var(z)
# Sum of the variances
3^2 + 4^2
# Compare to Normal distribution
curve(dnorm(x, mean = 30, sd = sqrt(3^2 + 4^2)), add = TRUE, col = "red", lwd = 2)


# Quotients of Normals is Cauchy

x <- rnorm(B, mean = 0, sd = 1)
y <- rnorm(B, mean = 0, sd = 1)
z <- x / y
par(mfrow = c(3, 1))
hist(x, probability = TRUE, xlim = c(-50, 50))
hist(y, probability = TRUE, xlim = c(-50, 50))
hist(z, probability = TRUE, xlim = c(-50, 50), breaks = 50, main = "Quotient of normals")
# Compare to Cauchy distribution
curve(dcauchy(x, location = 0, scale = 1), add = TRUE, col = "red", lwd = 2)
# Mean and variance of the quotient
mean(z)
var(z)
# Neither exists for Cauchy distribution
# But median and IQR do
median(z)
IQR(z)
# Median and IQR of Cauchy distribution
qcauchy(0.5, location = 0, scale = 1)
qcauchy(0.75, location = 0, scale = 1) -
  qcauchy(0.25, location = 0, scale = 1)
# Very close!
# Note: the histograms of the quotient will vary
# a lot because of extreme values