In [3]:

    

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np


    

In [53]:

    

from scipy.stats import norm, t

M = 20000 # number of monte carlo trials
N = 5 # number of samples, as if a household survey
mu = 0
sigma = 1

q50_err_nonpara = list()
q50_err_para = list()
q10_err_nonpara = list()
q10_err_para = list()
for m in range(M):
    #sample = norm.rvs(mu, sigma, size=N)
    sample = t.rvs(5, size=N)
    muhat = np.mean(sample)
    sigmahat = np.sqrt(np.var(sample))

    q50_err_nonpara.append(np.median(sample))
    q50_err_para.append(np.mean(sample))
    q10_err_nonpara.append(np.percentile(sample, 10))
    q10_err_para.append(norm.ppf(0.10,muhat,sigmahat))
    
print("True median:", t.ppf(0.5, 5))
print("True 10 percentile:", t.ppf(0.1, 5))


    

    




True median: 6.976003623e-17
True 10 percentile: -1.47588404878

In [54]:

    

plt.hist(q50_err_para, bins = 50, color="blue");
plt.hist(q50_err_nonpara, bins = 50, color="red");


    

In [55]:

    

plt.hist(q10_err_para, bins = 50, color="blue");
plt.hist(q10_err_nonpara, bins = 50, color="red");


    

In [ ]: