Babies 1: fig 1.8

# A function for calculating coefficient of kurtosis
kurt = function(x)
{
return(mean(((x-mean(x, na.rm=T))/sd(x, na.rm=T))^4))
}

isSmoker = bab1$smoke == 1
currKurt = kurt(bab1$bwt[isSmoker])

# simulate a normal distribution 1000 times
numSamples = sum(isSmoker, na.rm=T)
simDat = matrix(rnorm(1000*numSamples), ncol=1000)
simKurt = apply(simDat, 2, kurt)
hist(simKurt, xlab="Kurtosis value", ylab="Number of samples", xlim=c(2,4), main="")