oberon% splus
S-PLUS : Copyright (c) 1988, 1996 MathSoft, Inc.
S : Copyright AT&T.
Version 3.4 Release 1 for HP 9000 Series, HP-UX 9.0.5 : 1996 
Working data will be in /u/disk13/stat/mai/.Data 
> 
> 
> dim(result)
[1] 100   6
> ls()
  [1] ".Last.value"    ".Random.seed"   "BlindSimu"      "FFF"           
  [5] "KMNAatT"        "KMT"            "KMatT"          "KapMei"        
  [9] "LogRank2"       "Wdataclean"     "X"              "X1"            
 [13] "b"              "basefit"        "basefit1"       "basefit2"      
 [17] "bca1kexp25"     "bcanon"         "beta"           "bibj"          
 [21] "bibjEMreg"      "bj500exp2"      "bj500t2"        "bj500u"        
 [25] "bj500u2"        "boot2"          "c1"             "c2"            
 [29] "cancer"         "canstatus"      "cantime"        "cen"           
 [33] "coxfit"         "cxrg.cancr.fit" "d"              "d1"            
 [37] "dataclean"      "dd"             "del"            "delta"         
 [41] "finalbase"      "fq"             "fqout"          "fres"          
 [45] "fsortcum"       "g"              "gee"            "gp"            
 [49] "gtemp"          "gtemp.out"      "gwm"            "h"             
 [53] "hwtol"          "i1"             "i2"             "ind"           
 [57] "invreg"         "invreg.exp"     "invreg.h"       "invreg.log"    
 [61] "invreg.t"       "invreg.u"       "isot"           "iter"          
 [65] "jackknife"      "jump60"         "last.dump"      "last.warning"  
 [69] "lat"            "log500exp"      "logit500exp"    "loglik"        
 [73] "loklik"         "lsfittheta"     "mka"            "mres.ag.fit"   
 [77] "mu"             "myisot"         "newbj2000t"     "newt"          
 [81] "newtemp1"       "newtemp2"       "newtemp3"       "norm.inter"    
 [85] "num"            "perc"           "pred60"         "prim"          
 [89] "psum"           "result"         "rmp"            "rmp1"          
 [93] "roots"          "rx"             "ry"             "sd"            
 [97] "seedfor500"     "surv60"         "sx"             "temp1"         
[101] "temp12"         "temp503"        "temp5032"       "theta"         
[105] "time60"         "tsi"            "twosample"      "u"             
[109] "u500"           "unif250"        "x"              "x.cancer"      
[113] "x1"             "x12"            "x2"             "xbeta"         
[117] "xx"             "xxsort"         "xxx"            "y"             
[121] "y1"             "yy"             "yyy"            "z"             
> BlindSimu
function(X1 = rexp(50), C1 = rexp(50), X2 = rexp(50), C2 = rexp(50))
{
        Z1 <- pmin(X1, C1)
        D1 <- as.numeric(C1 > X1)
        Z2 <- pmin(X2, C2)
        D2 <- as.numeric(C2 > X2)
        temp1 <- KMNAatT(Z1, D1, 0.5)
        temp2 <- KMNAatT(Z2, D2, 0.5)
        Tvar <- (temp1$std.err)^2 + (temp2$std.err)^2
        temp3 <- KMNAatT(c(Z1, Z2), c(D1, D2), 0.5)
        Bvar <- 4 * (temp3$std.err)^2
        list(Tvariance = Tvar, Bvariance = Bvar, TnelsonVar = (temp1$
                NelsonAalenVar + temp2$NelsonAalenVar), BnelsonVar = 4 * (temp3$
                NelsonAalenVar), Best = temp3$surv, BGW = temp3$std.err)
}
> exp(-0.5)
[1] 0.6065307
> exp(-1)
[1] 0.3678794
> KMNAatT
function(z, d, t)
{
# This code actually computes a weighted Kaplan-Meier estimator at t.
# Z may have many ties. The function dataclean() sort the z and remove the
# tied observations and put the multiplicity in neww, the weight. 
# May be I should put a new input vector, weights, up there. 
# z can still have tie between a censored and uncensored obs.
# in this case the uncensored one should come first, i.e. in d, 1 comes first.
# d is of the same length of z.
# The last obs may need to be redefined to uncensored if you want a true
# distribution as estimator always.
# when the given t is larger then all z, then there may be a bias in the est.
# in this case the output will have a bias = T
        n <- length(z)
        if(n < 2)
                stop("times vector, z, is shorter then 2!")
        if(length(d) != n)
                stop("two inputs must be of same length")
        if(any((d != 0) & (d != 1)))
                stop("d must be 0/1's for censor/not-censor")
        temp <- dataclean(z, d)
        d <- temp$d
        n <- length(d)
        w <- temp$w
        sur <- rep(1, n)
        jum <- rep(0, n)
        var <- rep(0, n)
        nelson <- rep(0, n)
        nelsonVar <- rep(0, n)
        nelsonJum <- rep(0, n)
        if(d[1] == 1) {
                nelson[1] <- w[1]/sum(w)
                nelsonJum[1] <- nelson[1]
                nelsonVar[1] <- nelson[1]/(sum(w))
                var[1] <- w[1]/(sum(w) * (sum(w[2:n])))
                sur[1] <- 1 - nelson[1]
                jum[1] <- 1 - sur[1]
        }
        if(n < 2)
                stop("only one distinct time in z?")
        for(i in 2:n)
                if(d[i] == 0) {
                        var[i] <- var[i - 1]
                        sur[i] <- sur[i - 1]
                        nelson[i] <- nelson[i - 1]
                        nelsonVar[i] <- nelsonVar[i - 1]
                }
                else {
                        nelsonJum[i] <- w[i]/sum(w[i:n])
                        nelson[i] <- nelson[i - 1] + nelsonJum[i]
                        nelsonVar[i] <- nelsonVar[i - 1] + nelsonJum[i]/sum(w[i:
                                n])
                        var[i] <- var[i - 1] + nelsonJum[i]/(sum(w[(i + 1):n]))
                        sur[i] <- sur[i - 1] * (1 - nelsonJum[i])
                        jum[i] <- sur[i - 1] - sur[i]
                }
        sigma <- (sur) * sqrt(var)
        m <- sum(temp$z <= t)
        PosBias <- (t > max(z))
        if(PosBias)
                warning("the estimate is biased")
        list(time = t, surv = sur[m], std.err = sigma[m], NelsonAalen = nelson[
                m], NelsonAalenVar = nelsonVar[m], bias = PosBias)
}
> set.seed(153)
> set.seed(153)
> result <- matrix(NA, nrow=5000, ncol=6)
> for(i in 1:1000) result[i,] <- unlist(BlindSimu())
> for(i in 1001:2000) result[i,] <- unlist(BlindSimu())
> for(i in 2001:3000) result[i,] <- unlist(BlindSimu())
> for(i in 3001:4000) result[i,] <- unlist(BlindSimu())
> for(i in 4001:5000) result[i,] <- unlist(BlindSimu())
> 

> mean(result[,1]
+ )
[1] 0.01239244
> mean(result[,2])                   #Blind variance est
[1] 0.01252191
> mean(result[,1])                   #correct var est
[1] 0.01239244
> mean(result[,3])                   #correct nelson var est
[1] 0.03509194
> mean(result[,4])                   # blind nelson var est 
[1] 0.03476628
> var(result[,5])                    # sample variance , mixed sample
[1] 0.003142266
> mean(result[,6])    
[1] 0.05589574
> (mean(result[,6]))^2               # Greenwood, mixed sample
[1] 0.003124333
> 

or is this? 

> mean( (result[,6])^2 )
[1] 0.003130477