My aim is to study "Sample distributions" via simulations convincing myself of the outcomes.
After struggling with questions here, here and here, finally I hope am somewhat getting my head around it.
I tried to simulate a population of bernoulli random variable Y, with a given probability for Y=1. I tried to sample from it, each time sample size being 'n', for 'N' experiments. The sampling distribution turned out to be normal as expected.
I tried to then show, the problem would not be a good fit for normal approximation, when np < 10 or nq < 10 because the distribution would be skewed. But when I simulate in python, I do not observe much skewness. The "extreme" curves were as good as some of intermediate curves, so I do not get it.
I tried increasing sample size, no of experiments, but no improvement. I tried using bar instead of hist, but in vain. I tried varying width of them but no use (in fact changing width affects size for another purpose - to maintain total height as 1 due to density=True)
Kindly check and let me know if there is any issue in the math I am doing.
Code:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from SDSPSM import create_bernoulli_population, sample_for_SDSP
import numpy as np
from math import sqrt, pi
# control inputs
T = 10000
p_list = np.arange(0.1,1,0.1) # varying probability
N = 2000 # no of experiments
n = 50 # sample size
fig, ax = plt.subplots(1,1,figsize=(12,4))
plt.close()
def animate(i): ax.clear() p = round(p_list[i],4) # create population to be sampled from - Note pop has to be re created every time due to changing p pops = create_bernoulli_population(T, p) _, Y_mean_list = sample_for_SDSP(pops, N,n) # plot discrete density _, bins,_ = ax.hist(Y_mean_list, density=True) # normal approx mu, var, sigma = get_metrics(Y_mean_list) X = np.linspace(min(bins),max(bins),10*len(bins)) if sigma>0: Cp = 1/(sigma*sqrt(2*pi)) Ep = -1/2*((X-mu)/sigma)**2 G = Cp*np.exp(Ep) ax.plot(X, G, color='red') metrics_text = '$\mu_x:{}$ \n$\sigma_x:{}$'.format(mu, sigma) ax.text(0.97, 0.98,metrics_text,ha='right', va='top',transform = ax.transAxes,fontsize=10,color='red') ax.set_ylim([0,15]) ax.set_xlim([0,1]) np_factor = round(n*p,2) nq_factor = round(n*(1-p),2) metrics_text_2 = '$p:{}$ \n$np:{}$ \n$nq:{}$'.format(p, np_factor, nq_factor) ax.text(0.01, 0.98,metrics_text_2,ha='left', va='top',transform = ax.transAxes,fontsize=10,color='red')
plt.show()
ani1 = animation.FuncAnimation(fig, animate, frames = range(0,len(p_list)), interval=1000)
from IPython.display import HTML
def display_animation(anim): plt.close(anim._fig) return HTML(anim.to_html5_video())
display_animation(ani1)Helper file: SDSPSM
Output:
