set verbose off
clear

nulldata 200000
shape = 64  # gamma shape parameter
scale = 625 # gamma scale parameter
Income = randgen(G, shape, scale)
freq Income --normal --plot=display --silent
μ = mean(Income)
σ = sd(Income)
skew = skewness(Income)

printf "Parameters of the 'population' of %d incomes:\n", $nobs
printf " μ = %g, σ = %g, skewness %g\n\n", μ, σ, skew

K = 3000 # replications
n =  800 # sample size

printf "Theoretical σ(xbar) = σ/sqrt(n) = %g\n", σ/sqrt(n)

matrix xbar = zeros(K,1)
matrix ssd = zeros(K,1)

loop i=1..K
	smpl n --random --quiet
	xbar[i] = mean(Income)
	ssd[i] = sd(Income)
	smpl full --quiet
endloop

printf "\nStatistics based on the %d replications:\n", K
printf " Mean of xbar = %g\n", meanc(xbar)
printf " Mean of sample sd = %g\n", meanc(ssd)
printf " Mean of s_xbar = %g\n", meanc(ssd / sqrt(n))

printf "\nx-bar distribution:\n"
printf " sd %g, skewness %g\n", sd(xbar), skewness(xbar)

freq --matrix=xbar --normal --plot=display --silent