> #7.7 Radioactive Counts Example
> x<-rbind(
+ c(0,  57),
+ c(1, 203),
+ c(2, 383),
+ c(3, 525),
+ c(4, 532),
+ c(5, 408),
+ c(6, 273),
+ c(7, 139),
+ c(8,  45),
+ c(9,  27),
+ c(10, 10),
+ c(11,  4),
+ c(12,  0),
+ c(13,  1),
+ c(14,  1))
> 
> #Find mean to work out expected counts from Poisson distribution
> v<-sum(x[,1]*x[,2])/sum(x[,2])
> v
[1] 3.871549
> probs<-v^(x[,1])*exp(-v)/factorial(x[,1])
> probs[length(probs)]<-1-sum(probs[1:(length(probs)-1)])
> tab<-cbind(x[,2],sum(x[,2])*probs)
> colnames(tab)<-c("Observed","Expected")
> rownames(tab)<-0:(length(probs)-1)
> tab
   Observed    Expected
0        57  54.3144249
1       203 210.2809617
2       383 407.0565318
3       525 525.3131137
4       532 508.4438755
5       408 393.6930837
6       273 254.0336826
7       139 140.5005529
8        45  67.9943483
9        27  29.2492729
10       10  11.3239996
11        4   3.9855836
12        0   1.2858652
13        1   0.3829454
14        1   0.1417581
> 
> #totals & plot
> apply(tab,2,sum)
Observed Expected 
    2608     2608 
> 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)
> 
> #group so expected>1
> pooledtab<-tab
> while (any(pooledtab[,2]<=1))
+ {
+ n<-dim(pooledtab)[1]
+ pooledtab<-rbind(pooledtab[1:(n-2),],apply(pooledtab[(n-1):n,],2,sum))
+ }
> pooledtab
   Observed   Expected
0        57  54.314425
1       203 210.280962
2       383 407.056532
3       525 525.313114
4       532 508.443876
5       408 393.693084
6       273 254.033683
7       139 140.500553
8        45  67.994348
9        27  29.249273
10       10  11.324000
11        4   3.985584
          2   1.810569
> 
> #gof stuff - note have to specify how many parameters were estimated
> GOF(pooledtab,q=1)
Note 1 parameters were estimated to calculate expected values
The particular component indices given zero df to compensate are:
12 
This information is only used to calculate partition dfs and p-values

$AndersonDarling
[1] 1.248217

$Vr
 [1] -0.022771701 -1.806824765 -0.397543938  2.251815196  0.019319432
 [6] -0.710470278 -0.603115762  0.001398022  1.204697668  0.119051129
[11] -1.277152583 -0.717405930

$partition
           df           SS     pvalue
Location    1 5.185504e-04 0.98183238
Dispersion  1 3.264616e+00 0.07078961
Skewness    1 1.580412e-01 0.69096639
Kurtosis    1 5.070672e+00 0.02433395
Residual    7 4.480151e+00 0.72310604
Total      11 1.297400e+01 0.29502825

> 
> #gof stuff - attribute the estimated parameter to solving for the location component to be zero
> GOF(pooledtab,q=1,qcomps=c(1))
Note 1 parameters were estimated to calculate expected values
The particular component indices given zero df to compensate are:
1 
This information is only used to calculate partition dfs and p-values

$AndersonDarling
[1] 1.248217

$Vr
 [1] -0.022771701 -1.806824765 -0.397543938  2.251815196  0.019319432
 [6] -0.710470278 -0.603115762  0.001398022  1.204697668  0.119051129
[11] -1.277152583 -0.717405930

$partition
           df           SS     pvalue
Location    0 5.185504e-04         NA
Dispersion  1 3.264616e+00 0.07078961
Skewness    1 1.580412e-01 0.69096639
Kurtosis    1 5.070672e+00 0.02433395
Residual    8 4.480151e+00 0.81141559
Total      11 1.297400e+01 0.29502825

>