Nagarwalla.test<-function(time,t,n0,Nrep){
# ---------------------------------------------------------------
# An example of execution
#
#    yday<-scan("d:/demo/Knox.dat")
#    n0<-5; tt<-2191; Nrep<-9999
#    out<-Nagarwalla.test(yday,tt,n0,Nrep)
#
# ---------------------------------------------------------------
# ARGUMENTS
#
#   time : a vector of individual time points
#   t    : T, length of study period
#   n0   : Nagarwalla's predetermined value of n0
#   Nrep : Number of Monte Calro replications, e.g., 999, 9999
#
# VALUES
#
#   Nagarwalla.T       : test statistic
#   kstar              : the number k of the points in the cluster
#   MLCfrom            : starting time points of the cluster
#   MLCto              : terminal time points of the cluster
#   Freq               : a vector of simulated test statistics under the null
#   Simulated.p.value  : simulated p-value
# ---------------------------------------------------------------
xday<-time / t
nt<-length(time)
#
dmin<-rep(0,nt)
Imin<-dmin
lg<-dmin
#
for (k in 1:(nt-1-n0) ){
   n<- n0-1+k
   x1<-1:(nt-n+1)
   for (i in 1:(nt-n+1) ){
      x1[i]<- xday[i+n-1]-xday[i]
}
   dmin[n]<-min(x1)
   dxx    <-min(x1)
   ind<- 1:length(x1)
#   Imin[n]<-ind[rank(x1)==1]
  for (i in 1:(nt-n+1) ){
    if (abs(x1[i]-dxx)<0.00001){ jj<-i }
}
   Imin[n]<-jj
   lg[n]<- n*log(n/nt) + (nt-n)*log((nt-n)/nt) - 
       n*log(dmin[n]) - (nt-n)*log(1-dmin[n])
}
   inds<- 1:length(lg)
   ns<-inds[rank(lg)==length(lg)]
#   ans<-c(ns, Imin[ns], exp(lg[ns]))
   jjns<-ns
   jjj1<-Imin[ns]
   jjj2<-jjj1+ns-1
obst<-lg[ns]
#
# Monte Carlo Hypothesis testing 
#
ktmon<-1:(Nrep+1)
for(j in 1:Nrep){
xday<-sort(runif(nt))
#
dmin<-rep(0,nt)
Imin<-dmin
lg<-dmin
#
for (k in 1:(nt-1-n0) ){
   n<- n0-1+k
   x1<-1:(nt-n+1)
   for (i in 1:(nt-n+1) ){
      x1[i]<- xday[i+n-1]-xday[i]
}
   dmin[n]<-min(x1)
   dxx    <-min(x1)
   ind<- 1:length(x1)
#   Imin[n]<-ind[rank(x1)==1]
  for (i in 1:(nt-n+1) ){
    if (abs(x1[i]-dxx)<0.00001){ jj<-i }
}
   Imin[n]<-jj
   lg[n]<- n*log(n/nt) + (nt-n)*log((nt-n)/nt) - 
       n*log(dmin[n]) - (nt-n)*log(1-dmin[n])
}
   inds<- 1:length(lg)
   ns<-inds[rank(lg)==length(lg)]
#
ktmon[j] <- lg[ns]
}
ktmon[Nrep+1]<-obst
r<-length(ktmon[ktmon>=obst])
p<-r/(Nrep+1)
list(Nagarwalla.T=obst, kstar=jjns, MLCfrom=time[jjj1], MLCto=time[jjj2], 
Freq=ktmon, Simulated.p.value=p)
}