#problem 45 chapter 9
#number of male children in 6115 families with 12 children
z = read.csv("geissler.txt")
number = z[,1]
freq = z[,2]  #observed counts
p1 = dbinom(0:12,12,.5)  #binomial probs
exp = 6115*p1 #expected counts
plot(0:12,freq)
points(0:12,exp,pch=19)
plot(0:12, freq-exp)
x2 = sum((freq-exp)^2/exp)
tst = chisq.test(freq,p=p1)
tst$statistic
1-pchisq(tst$statistic, 12)  #pvalue

#now test with phat estimated by max likelihood
phat = sum(number*freq/(6115*12)) #number of males/total number children
p2 = dbinom(0:12,12,phat)
exp = 6115*dbinom(0:12,12,phat)
plot(0:12,freq - exp)
x2 =chisq.test(freq, p=p2)  #df is wrong,should be 11
x2$statistic
1 - pchisq(x2$statistic,11)  #p-value
plot(0:12, x2$residuals) # residuals = (obs -exp)/sqrt(expected)