set.seed(1) N=100 #X=rexp(N) #Y=rnorm(N) #qqnorm(Y) #qqnorm(X) b=exp(6) g=1/3 y=rweibull(100,g,b) c=runif(100,0,1780) d=as.numeric(y<=c) m=y*d+c*(1-d) #in R library(splines) library(survival) #survreg in R =survReg in Splus #x<-matrix(scan("junk"), ncol=2, byrow=T) x=sort(m) #u=survfit(Surv(m,d)~1) #plot(u,lty=1) (zz=survreg(Surv(m,d)~1,dist="gaussian")) plot(x,1-pnorm(x,zz$coef,zz$scale),type="l",lty=2,xlim=c(0,1500)) (zz=survreg(Surv(m,d)~1,dist="loggaussian")) lines(x,1-plnorm(x,zz$coef,zz$scale),type="l",lty=5) (zz=survreg(Surv(m,d)~1)) lines(x,1-pweibull(x,1/zz$scal,exp(zz$coef)),type="l",lty=3) r=mean(d)/mean(m) lines(x,exp(-r*x),type="l",lty=4) leg.names<-c("norm","weib","exp","lognormal") legend(300, 0.78, leg.names, lty=c(2,3,4,5),cex=1.0) u=survfit(Surv(m,d)~1) plot(u,lty=1,xlim=c(0,1500)) x=sort(m) (zz=survreg(Surv(m,d)~1,dist="gaussian")) lines(x,1-pnorm(x,zz$coef,zz$scale),type="l",lty=2) (zz=survreg(Surv(m,d)~1,dist="lognormal")) lines(x,1-plnorm(x,zz$coef,zz$scale),type="l",lty=5) (zz=survreg(Surv(m,d)~1)) lines(x,1-pweibull(x,1/zz$scal,exp(zz$coef)),type="l",lty=3) r=mean(d)/mean(m) lines(x,exp(-r*x),type="l",lty=4) leg.names<-c("ple", "norm","weib","exp","lognormal") legend(300, 0.88, leg.names, lty=c(1,2,3,4,5),cex=1.0) u=survfit(Surv(m,d)~1) #plot(u,lty=1) plot(u$time, u$surv, type="s",lty=1, xlim=c(0,1500)) x=sort(m) #(zz=survreg(Surv(m,d)~1,dist="gaussian")) #lines(x,1-pnorm(x,zz$coef,zz$scale),type="l",lty=2) (zz=survreg(Surv(m,d)~1,dist="lognormal")) (Z=summary(zz)) Z$table (w=Z$table[1,2]) #lines(x,1-plnorm(x,zz$coef,zz$scale),type="l",lty=5) lines(x,1-plnorm(x,zz$coef-1.96*w,zz$scale),type="l",lty=5) lines(x,1-plnorm(x,zz$coef+1.96*w,zz$scale),type="l",lty=5) (zz=survreg(Surv(m,d)~1)) #lines(x,1-pweibull(x,1/zz$scal,exp(zz$coef)),type="l",lty=3) (Z=summary(zz)) Z$table (w=Z$table[1,2]) lines(x,1-pweibull(x,1/zz$scal,exp(zz$coef-1.96*w)),type="l",lty=3) lines(x,1-pweibull(x,1/zz$scal,exp(zz$coef+1.96*w)),type="l",lty=3) #r=mean(d)/mean(m) #lines(x,exp(-r*x),type="l",lty=4) leg.names<-c("ple", "weib","lognormal") legend(300, 0.88, leg.names, lty=c(1,3,5),cex=1.0) q() set.seed(1) N=100 X=rexp(N) Y=rnorm(N) qqnorm(Y) qqnorm(X) b=exp(6) g=1/3 y=rweibull(100,g,b) c=runif(100,0,780) d=as.numeric(y<=c) m=y*d+c*(1-d) #in R library(splines) library(survival) #survreg in R =survReg in Splus #x<-matrix(scan("junk"), ncol=2, byrow=T) m d #in R library(splines) library(survival) #survreg in R =survReg in Splus #x<-matrix(scan("junk"), ncol=2, byrow=T) x=sort(m) #u=survfit(Surv(m,d)~1) #plot(u,lty=1) (zz=survreg(Surv(m,d)~1,dist="gaussian")) plot(x,1-pnorm(x,zz$coef,zz$scale),type="l",lty=2) #(zz=survreg(Surv(log(m),d)~1,dist="gaussian")) #lines(x,1-pnorm(log(x),zz$coef,zz$scale),type="l",lty=5) (zz=survreg(Surv(m,d)~1,dist="loggaussian")) lines(x,1-plnorm(x,zz$coef,zz$scale),type="l",lty=5) lines(x,1-pnorm(log(x),zz$coef,zz$scale),type="l",lty=1) (zz=survreg(Surv(m,d)~1)) lines(x,1-pweibull(x,1/zz$scal,exp(zz$coef)),type="l",lty=3) #lines(x,1-pnorm(log(x),zz$coef,zz$scale),type="l",lty=5) #it used incorrectly weib parameters. #(zz=survreg(Surv(m,d)~1,dist="exponential")) r=mean(d)/mean(m) lines(x,exp(-r*x),type="l",lty=4) #u=survfit(Surv(m,d)~1) #lines(u,lty=1) leg.names<-c("norm","weib","exp","lognormal") legend(300, 0.78, leg.names, lty=c(2,3,4,5),cex=1.0) u=survfit(Surv(m,d)~1) plot(u,lty=1) x=sort(m) (zz=survreg(Surv(m,d)~1,dist="gaussian")) lines(x,1-pnorm(x,zz$coef,zz$scale),type="l",lty=2) #lines(x,1-pnorm(log(x),zz$coef,zz$scale),type="l",lty=1) (zz=survreg(Surv(m,d)~1,dist="loggaussian")) #lines(x,1-plnorm(x,zz$coef,zz$scale),type="l",lty=5) (zz=survreg(Surv(log(m),d)~1,dist="gaussian")) lines(x,1-plnorm(x,zz$coef,zz$scale),type="l",lty=5) (zz=survreg(Surv(m,d)~1)) lines(x,1-pweibull(x,1/zz$scal,exp(zz$coef)),type="l",lty=3) #(zz=survreg(Surv(m,d)~1,dist="exponential")) r=mean(d)/mean(m) lines(x,exp(-r*x),type="l",lty=4) leg.names<-c("ple", "norm","weib","exp","lognormal") legend(300, 0.88, leg.names, lty=c(1,2,3,4,5),cex=1.0) (zz=survreg(Surv(m,d)~1,dist="lognormal")) set.seed(1) x=rexp(100,2) y=runif(100,0,4) d=as.numeric(x<=y) z=x*d+y*(1-d) u=survfit(Surv(m,d)~1,conf.int=0.99) plot(u$time,u$surv,type="s",lty=1) x=sort(m) plot(x,exp(-2*x),type="l",lty=1) (zz=survreg(Surv(y)~1, dist="gaussian")) z=rexp(100) y=exp(6+3*log(z)) # log y=6+3*ln z mean(y) #(zz=survReg(Surv(y)~1, dist="exponential")) #(zz=survReg(Surv(y)~1)) #(zz=survReg(Surv(y)~1, dist="weibull")) sqrt(pi/100) y=rweibull(100,2, 1/5) mean(y) set.seed(1) y=rweibull(100,1/3, exp(6)) (zz=survreg(Surv(y)~1)) summary(zz) set.seed(1) b=exp(6) g=1/3 y=rweibull(100,g,b) c=runif(100,0,780) d=as.numeric(y<=c) m=y*d+c*(1-d) (zz=survreg(Surv(m,d)~1)) #summary(zz) #plot(zz) #summary(aa) #plot(aa) #Monte Carlo set.seed(1) N=200 L=0 M=c(0,0) for(i in 1:N) { bb=rnorm(1,b,100) gg=rnorm(1,g,1/9) a=prod(d*dweibull(m,gg,bb)+(1-d)*(1-pweibull(m,gg,bb))) if (L