################################################################################

# CLT

# Normal

n <- 30		# sample size
k <- 1000		# number of samples

mu <- 5		# population mean, mu, for the "unknown" population
sigma <- 2		# population standard deviation, sigma, for the "unkown" population

x.down <- mu - 3*sigma				# gives the lower x-value for the population
x.up <- mu + 3*sigma					# gives the upper x-value for the population

x <- seq(x.down,x.up,0.01)	

y.up <- 1.5 	# give the maximum value for the y-axis

y <- dnorm(x,mu,sigma)		# calculates the normal distribution.

plot(x,y,type='l',xlim=c(x.down,x.up),ylim=c(0,y.up), main='Normal population')

SEM <- sigma/sqrt(n)		# standard error of the mean

x <- matrix(rnorm(n*k,mu,sigma),n,k)		# This gives a matrix with k = 1000 samples down the columns.

x.mean <- apply(x,2,mean)		# computes the means down the columns for the k samples

hist(x.mean,prob=T,xlim=c(x.down,x.up),ylim=c(0,y.up),main='Sampling distribution of the sample mean, Normal case')

par(new=T)			# resets the picture so another picture can be plotted on the same graph.

x <- seq(x.down,x.up,0.01)		

y <- dnorm(x,mu,SEM)		# calculates the theoretical normal distribution.

plot(x,y,type='l',xlim=c(x.down,x.up),ylim=c(0,y.up))

# Exponential

n <- 30		# sample size
k <- 1000		# number of samples

theta <- 3		# mean of the Exponential

lambda <- 1/theta

mu <- theta
sigma <- theta

x.down <- 0.01 #mu - 3*sigma				# gives the lower x-value for the population
x.up <- mu + 3*sigma					# gives the upper x-value for the population

x <- seq(x.down,x.up,0.01)		

y.up <- 1.0

y <- dexp(x,lambda)		# calculates the Exponential distribution.

plot(x,y,type='l',xlim=c(x.down,x.up),ylim=c(0,y.up), main='Exponential population')

SEM <- sigma/sqrt(n)

x <- matrix(rexp(n*k,lambda),n,k)		# This gives a matrix with the samples down the columns.

x.mean <- apply(x,2,mean)

hist(x.mean,prob=T,xlim=c(x.down,x.up),ylim=c(0,y.up),main='Sampling distribution of the sample mean, Normal case')

par(new=T)

x <- seq(x.down,x.up,0.01)

y <- dnorm(x,mu,SEM)

plot(x,y,type='l',xlim=c(x.down,x.up),ylim=c(0,y.up))