Stone.test<-function(oo, ee, Nrep){
# -----------------------------------------------------------------
#  An example of execution
#
#     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<-Stone.test(oo,ee,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
#     Nrep    : The number of Monte Carlo replications, e.g., 999, 9999
#
#  VALUES
#     Test.statistic : test statistic
#     Estimate       : a vector of maximum likelihood estimates of relative risk
#     p.value        : simulated p-value
#     Freq           : a vector of simulated test statistics under the null
#
#-----------------------------------------------------------------------
nr<-length(oo)
gg<-1:nr; g<-gg; xlam<-gg
#
for (i in 1:nr){
  gg<-rep(10000000,nr)
  for (is in 1:i){
    g<-rep(0,nr) 
    for (it in i:nr){
      os<-0; es<-0
      for (ir in is:it){
        os<-os + oo[ir]
        es<-es + ee[ir]
      }
      g[it]<-os/es
    }
    smax<- -20
    for (j in i:nr){
      if (g[j] >= smax) { smax<- g[j] }
    }
    gg[is]<- smax
  }
  smin<-10000
  for (j in 1:i){
    if (gg[j] <= smin) { smin<- gg[j] }
  }
  xlam[i]<- smin
}
mean<- sum(xlam*ee)/sum(ee)
Ts <- sum( oo * (log(xlam) - log(mean)) )
Est<- xlam
plot(1:nr, xlam)
#
nn<-sum(oo)
pvec<-ee/sum(ee)
nloop<-Nrep+1
for (i in 2:nloop){
aa<-rmultinom(1,nn,pvec)
#
for (i in 1:nr){
  gg<-rep(10000000,nr)
  for (is in 1:i){
    g<-rep(0,nr) 
    for (it in i:nr){
      os<-0; es<-0
      for (ir in is:it){
        os<-os + aa[ir]
        es<-es + ee[ir]
      }
      g[it]<-os/es
    }
    smax<- -20
    for (j in i:nr){
      if (g[j] >= smax) { smax<- g[j] }
    }
    gg[is]<- smax
  }
  smin<-10000
  for (j in 1:i){
    if (gg[j] <= smin) { smin<- gg[j] }
  }
  xlam[i]<- smin
}
mean<- sum(xlam*ee)/sum(ee)
Ts1 <- sum( aa * (log(xlam) - log(mean)) )
Ts<- c(Ts, Ts1)
}
p1<- rank(-Ts)[1]/(Nrep+1)
list(Test.statistic=Ts[1], Estimate=Est, p.value=p1, Freq=Ts)
}