notebook.community

Edit and run



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)) 
lines(density(s.10, adjust=2), lty="dotted")  
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)) 
lines(density(s.100, adjust=2), lty="dotted") 
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)) 
lines(density(s.1000, adjust=2), lty="dotted") 
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)) 
lines(density(s.2000, adjust=2), lty="dotted") 
mean(s.2000)
var(s.2000)









    Out[60]:




1.96405148426902






    Out[60]:




0.199093917595788



In [ ]: