library(MASS)
library(fitdistrplus)

k <- n <- x <- y <- z <- i <- s <- h <- NULL

#Simulation condition
k <- 1000 #Number of segment
n <- 1000 #Number of sample
logmean <- log(100) #fixed scale parameter of lognormal
logsigma <- 1.0 # shape parameter of lognormal

s <- numeric(n)

#Random variable generation (Lognormal)
for (i in 1:n) {
　 t <- rlnorm(k,logmean,logsigma)
　s[i] <- min(t)
}

#Random variable generation (Weibull)
#for (i in 1:n) {
#　t <- rweibull(k,beta0,eta0)
#　s[i] <- min(t)
#}

i <- NULL

x <- c(s[])
y <- sort(x)
xmin <- min(y)/5
xmax <- max(y)*5

cat("Number of data \n")
print(n)
cat(" \n")

#Non-parametric estimation of cumulative failure
for(i in 1:n){
　F <- (i-0.3)/(n+0.4) #Meddian rank
　#F <- i/(n+1) #Mean rank
　z <- c(z,F)
}

#Grid in figure
y0 <- c(10^(-10:0),0.1,0.5,1,5,10,20,30,40,50,60,70,80,90,95,99,99.9)
probs <- log(log(1/(1-(y0/100))))

#Preprocessing for plot
Y1 <- log(log(1/(1-z)))
X1 <- log(y)
Ymin <- min(Y1)-1
Ymax <- max(Y1)+0.5

#Parameter estimation by Weibull plot
result2 <- lsfit(Y1, X1, wt=NULL, intercept=TRUE, tolerance=1e-7, yname=NULL)
cf <- as.numeric(coef(result2))
beta <- 1/cf[2]
eta <- exp(cf[1])
corr <- cor(X1,Y1)

#Plotting
par(family="serif")
plot(y,Y1,log="x",xlim=c(xmin,xmax),ylim=c(Ymin,Ymax),xlab="time",ylab="Cumulative probability [%]",yaxt="n")

#Line fitting
lx <- c(xmin,xmax)
Yestmin <- beta*log(xmin)-beta*log(eta)
Yestmax <- beta*log(xmax)-beta*log(eta)
ly <- c(Yestmin,Yestmax)
lines(lx,ly)

#Grid
abline(h=probs,col="gray",lwd=1,lty=3)
axis(2,at=probs,labels=y0)
para <- c(beta, eta, corr)

#Parameter estimators by Weibull plot
cat("Estimators by Weibull plot \n")
cat(" Shape Scale C.C. \n")
print(para)
cat(" \n")

#Maximum likelihood method
wfit <- fitdist(y, "weibull")
result1 <- summary(wfit)
lnfit <- fitdist(y, "lnorm")
result2 <- summary(lnfit)
print(result1)
cat(" \n")
print(result2)

#Goodness of fit
cat(" \n")
result3 <- gofstat(list(wfit, lnfit), fitnames=c("Weibull", "lognormal"))
print(result3)