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<a)
{
L=a
M=c(gg,bb)
}
}
a
L
c(1/M[1],log(M[2]))
(x= c(zz$scale, zz$coef))
gg=1/x[1]
bb=exp(x[2])
prod(d*dweibull(m,gg,bb)+(1-d)*(1-pweibull(m,gg,bb)))

z=summary(zz)
names(z)
z$table
list(z)

# MLE of U(\theta,5),  with Y~U(0,4)
set.seed(1)
n=100
m=runif(n,0,5)
c=runif(n,0,4)
d=as.numeric(m<=c)
m=m*d+(1-d)*c
k=100
v=1:k
(t=min(m[d==1]))
for (i in 1:k){
m=runif(n,t,5)
c=runif(n,0,4)
d=as.numeric(m<=c)
m=m*d+(1-d)*c
m
d
(v[i]=min(m[d==1]))
}
sd(v)
t
v=sort(v)
v[5]
v[95]
v[1]
q()
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
     h8 <- hilbert(8); h8
     sh8 <- solve(h8)
     round(sh8 %*% h8, 3)
     
     A <- hilbert(4)
     A[] <- as.complex(A)
     ## might not be supported on all platforms
     try(solve(A))

q()