Bithell.test<-function(oo, ee, g, Nrep){ # ----------------------------------------------------------------- # An example of execution # # d<-1:10; g<- 1/d; Nrep<-999 # ob<- c(0, 1, 3, 2, 1, 8, 7, 4, 0, 2) # ex<- c(0.041,0.701,1.128,1.634,1.868,5.738,4.550,3.617,3.792,4.931) # out<-Bithell.test(oo,ee,g,Nrep) # # --------------------------------------------------------------------- # ARGUEMENTS # oo : a vector of observed number of cases sorted in the ascending # order of the distance to the point source # ee : a vector of conditional expected number of cases sorted # in the ascending order of the distance to the point source # g : a vector of values of exposure function of distance to # the point source # Nrep : The number of Monte Carlo replications, e.g., 999, 9999 # # VALUES # Test.statistic : test statistic # p.value : simulated p-value # Freq : a vector of simulated test statistics under the null # #----------------------------------------------------------------------- nr<-length(oo) Tb<- sum( oo * log(g) ) # nn<-sum(oo) pvec<-ee/sum(ee) nloop<-Nrep+1 for (i in 2:nloop){ aa<-rmultinom(1,nn,pvec) # Tb1<- sum( aa * log(g) ) Tb<- c(Tb, Tb1) } p1<- rank(-Tb)[1]/(Nrep+1) list(Test.statistic=Tb[1], p.value=p1, Freq=Tb) }