library(splines)
library(survival)
#4.7.2 Problem 2

rm(list=ls())
p = 1/2
n = 100
gmle1=rep(0,1000)
for(i in 1:1000){
x = rexp(n,p)
y1 = rexp(n,p)
z = rexp(n,p)
y2 = y1+z
status = 2*as.numeric(x<=y1)+3*as.numeric(y1<x)*as.numeric(x<=y2)+0*as.numeric(y2<x)
t1 = 0*(status==2)+y1*(status==3)+y2*(status==0)
t2 = y1*(status==2)+y2*(status==3)+100*(status==0)
my.surv = Surv(time=t1,time2=t2,event=status,type="interval")
temp = survfit(my.surv~1)
tt=summary(temp)
gmle1[i]=tt$surv[tt$time>=2][1]
}
s100=sd(gmle1)

p = 1/2
n = 400
gmle2=rep(0,1000)
for(i in 1:1000){
x = rexp(n,p)
y1 = rexp(n,p)
z = rexp(n,p)
y2 = y1+z
status = 2*as.numeric(x<=y1)+3*as.numeric(y1<x)*as.numeric(x<=y2)+0*as.numeric(y2<x)
t1 = 0*(status==2)+y1*(status==3)+y2*(status==0)
t2 = y1*(status==2)+y2*(status==3)+100*(status==0)
my.surv = Surv(time=t1,time2=t2,event=status,type="interval")
temp = survfit(my.surv~1)
tt=summary(temp)
gmle2[i]=tt$surv[tt$time>=2][1]
}
s400=sd(gmle2)
(d1=s100/s400-2)
(d2=s100/s400-4^(1/3))
(d3=s100/s400-4^(1/3)*(log(400)/log(100))^(1/3))



#4.7.2 Problem 3
rm(list=ls())
p = 1/5
n = 100
gmle1=rep(0,1000)
for(i in 1:1000){
x = rexp(n,p)
y1 = sample(1:9,n,replace=TRUE)
y2 = y1+1
status = 2*as.numeric(x<=y1)+3*as.numeric(y1<x)*as.numeric(x<=y2)+0*as.numeric(y2<x)
t1 = 0*(status==2)+y1*(status==3)+y2*(status==0)
t2 = y1*(status==2)+y2*(status==3)+1000*(status==0)
my.surv = Surv(time=t1,time2=t2,event=status,type="interval")
temp = survfit(my.surv~1)
tt=summary(temp)
gmle1[i]=tt$surv[tt$time>=2][1]
if(is.na(gmle1[i])) gmle1[i]=0
}
s100=sd(gmle1)

p = 1/5
n = 400
gmle2=rep(0,1000)
for(i in 1:1000){
x = rexp(n,p)
y1 = sample(1:9,n,replace=TRUE)
y2 = y1+1
status = 2*as.numeric(x<=y1)+3*as.numeric(y1<x)*as.numeric(x<=y2)+0*as.numeric(y2<x)
t1 = 0*(status==2)+y1*(status==3)+y2*(status==0)
t2 = y1*(status==2)+y2*(status==3)+100*(status==0)
my.surv = Surv(time=t1,time2=t2,event=status,type="interval")
temp = survfit(my.surv~1)
tt=summary(temp)
gmle2[i]=tt$surv[tt$time>=2][1]
if(is.na(gmle1[i])) gmle2[i]=0
}
s400=sd(gmle2)
(d1=s100/s400-2)
(d2=s100/s400-4^(1/3))
(d3=s100/s400-4^(1/3)*(log(400)/log(100))^(1/3))