Overview

In this project I will simulate the exponential distribution and compare my results to the Central Limit Theorem.

Simulations

#set up the simulation
n = 40 #length of exponetial
lambda = 0.2
mean = 1/lambda
sd = 1/lambda
nsim = 1000 # number of simulations 

#simulate e and convert e to a matrix
set.seed(2016)
e <- rexp(n*nsim,lambda)
e = matrix(e, nrow = nsim, ncol = n)

#calculate the averages of 1000 simulations of 40 random exponentials
emeans <- apply(e,1,mean)

Sample mean versus Theoretical mean

#calculate sample mean as the expected value of the distribution of 1000 averages
sample_mean <- mean(emeans)

#plot sample mean versus theoretical mean
barplot(c(sample_mean, mean), names.arg=c("sample mean", "theoretical mean"), col = c("red", "blue"))
title(main="Sample mean vs. theoretical mean of exponential distribution")
legend("topright", legend = c(round(sample_mean,2) , mean), col = c("red", "blue"), lty=c(1,1))

The sample mean of the expoential distribution is 4.98 while its theoretical mean is 5.

Sample variance versus Theoretical variance

#calculate the variance of 1000 simulations of 40 random exponentials
evars <- apply(e, 1, var)

#sample variance is the expected value of the distribution of 1000 variances 
sample_var <- mean(evars) 

#plot sample variance versus theoretical variance
barplot(c(sample_var, sd^2), names.arg=c("sample variance", "theoretical variance"), col=c("red", "blue"))
title(main="Sample variance vs. theoretical variance of exponential distribution")
legend("topright", legend = c(round(sample_var,2) , sd^2), col = c("red", "blue"), lty=c(1,1))

The sample variance of the expoential distribution is 24.93 while its theoretical mean is 25.

Distirbution (Central Limit Theorem)

#show the distribution of averages of 40 random exponentials is approximately normal
#first plot
#show the distribution of 40 random exponetials
par(mfrow = c(1,2), mar = c(4,2,4,2), cex.main = 1)
hist(rexp(40,lambda), main = "", xlab = "")
title(main = "Distribution of 40 random exponentials")

#second plot
#show the distribution of 1000 averages of 40 random exponentials
hist(emeans, main = "", xlab = "")
title(main = "Distribution of 1000 averages of 40 random exponentials", xlab="")

dev.off()
## null device 
##           1

The distribution of 1000 averages of 40 random exponentials looks more Gaussian or normal than the distribution of 40 random exponentials.