In [19]:

set.seed(124)




In [32]:

## A function to randomly sample 'sampling.size' points from
## 'sampling.distribution' 'replicates' number of times
## and store the means.

sampler <- function(replicates, sampling.size, sampling.distribtution){
sample.means <- vector()
N <- length(sampling.distribtution)
for(i in 1:replicates)
{
s <- sample(seq(1, N), sampling.size)
sample <- sampling.distribtution[s]
sample.means <- c(sample.means, mean(sample))

}
return(sample.means)
}




In [61]:

x.gamma <-rgamma(n=5000,shape=1,scale=2)
hist(x.gamma)
mean(x.gamma)
var(x.gamma)




Out[61]:

1.98301435580402

Out[61]:

4.06623608977553



For $x \sim \Gamma(a,b)$:

$$EX = ab$$

Hence the mean in our case should be 2.

$$Var = ab^2=4$$

And hence the sampling distribution of the mean should have $$Var(\bar{X}) = ab^2/n$$



In [56]:

s.10 <- sampler(10,20,x.gamma)
hist(s.10, main="Sampling distribution of the mean with 10 samples", xlab = "Sample mean", ylab="Probability", prob=T)
lines(density(s.10))
mean(s.10)
var(s.10)




Out[56]:

1.9144273721788

Out[56]:

0.370371464421344




In [57]:

s.100 <- sampler(100,20,x.gamma)
hist(s.100, main="Sampling distribution of the mean with 100 samples", xlab = "Sample mean", ylab="Probability", prob=T)
lines(density(s.100))
mean(s.100)
var(s.100)




Out[57]:

1.91777240446569

Out[57]:

0.154173668232302




In [58]:

s.1000 <- sampler(1000,20,x.gamma)
hist(s.1000, main="Sampling distribution of the mean with 1000 samples", xlab = "Sample mean", ylab="Probability", prob=T)
lines(density(s.1000))
mean(s.1000)
var(s.1000)




Out[58]:

1.94581638392367

Out[58]:

0.202825893452269




In [60]:

s.2000 <- sampler(2000,20,x.gamma)
hist(s.2000, main="Sampling distribution of the mean with 2000 samples", xlab = "Sample mean", ylab="Probability", prob=T)
lines(density(s.2000))
mean(s.2000)
var(s.2000)




Out[60]:

1.96405148426902

Out[60]:

0.199093917595788




In [ ]: