Overview

The goal of this exercise is to study exponential distribution and compare the result with Central Limit Theorem (CLT). The distribution of 1000 simulation of averages of 40 exponentials, and the theoretical distribution derived from CLT will be compared.

Simulations

Simulate 1000 times of averages of 40 exponentials, which is stored in vector mexp for following analysis.

simulations <- 1000
n <- 40
lambda <- 0.2
set.seed(7654)
x <- rexp(simulations, lambda)
mexp <- NULL
for (i in 1 : simulations) mexp <- c(mexp, mean(rexp(n, lambda)))

Sample Mean versus Theoretical Mean

As the lambda rate is 0.2 in this case, the theoretical mean is \(1/\lambda=5\). Calculate the sample mean.

mean(mexp)
## [1] 5.030255

It seems that the sample mean is quite close to the theoretical mean. Make a histogram to show the relation.

hist(mexp,main = "Sample Mean vs Theoretical Mean")
abline(v=mean(mexp),col="blue")
abline(v=5,col="red")
text(4, 240, paste("Sample mean = ", round(mean(mexp),2) ), col="blue")
text(5.7, 240, paste("Theoretical mean = 5" ), col="red")

Sample Variance versus Theoretical Variance

The theoretical standard deviation from CLT is \((1/\lambda)/\sqrt{n}\). Compare it with the sample standard deviation. The the variance is the square of them respectively.

print (paste("Theoretical variance: ", (1/lambda/sqrt(n))^2, ", Sample variance", round(var(mexp),3)))
## [1] "Theoretical variance:  0.625 , Sample variance 0.581"

Distribution

Plot the probability density of sample, and also plot the theoretical normal distribution with the theoretical mean and standard deviation. Compare two distribution, the sample distribution is approximately the theoretical distribution.

xnorm <- seq(3,9,length=1000)
ynorm <- dnorm(xnorm, mean=1/lambda, sd=1/lambda/sqrt(n))
hist(mexp, prob=TRUE, main="Distribution of averages of exponentials", xlab="Averages", breaks = 20)
lines(density(mexp), col="blue", lwd=2)
lines(xnorm, ynorm, type="l", col="grey", lty=2, lwd=3)
legend("topright", c("Sample distribution","Theoretical distribution"), lty=c(1,2), lwd=c(2,3),col=c("blue","grey"))