#7.8 Milk Bacteria Example
x<-rbind(
c(0,  56),
c(1, 104),
c(2,  80),
c(3,  62),
c(4,  42),
c(5,  27),
c(6,   9),
c(7,   9),
c(8,   5),
c(9,   3),
c(10,  2),
c(11,  1))

#Negative-binomial class probabilities from Jarvis (1989, p50)
#note the table on page 7.2 is incorrect
expected<-c(60.88,90.76,85.16,64.12,42.32,25.58,14.49,7.85,4.09,2.07,1.02,1.66)
tab<-cbind(x[,2],expected)
colnames(tab)<-c("Observed","Expected")
rownames(tab)<-x[,1]
tab

#totals & plot
apply(tab,2,sum)
barplot(t(tab),beside=T,legend.text=c("Observed","Expected"))

#DDR and plot
ddr<-sqrt(4*tab[,1]+2)-sqrt(4*tab[,2]+2)
ddr[tab[,1]==0]<-1-sqrt(4*tab[,2]+1)[tab[,1]==0]
plot(x[,1],ddr,xlab="Category",ylab="DDR",type="n")
points(x[,1],ddr,cex=2)

#gof stuff - note have to specify how many parameters were estimated
#Note here attribute the estimated parameters to solving for the location/scale components to be zero
GOF(tab,q=2,qcomps=c(1,2))