library(MASS)
library(fitdistrplus)

weiplot<-function(s){

x <- NULL
y <- NULL
z <- NULL
i <- NULL

x <- read.csv("weidata.csv", header=T)
x <- c(x[,s])
y <- sort(x)
xmin <- min(y)/5
xmax <- max(y)*5

n <- length(y)
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))))

#Linearity check
Y1 <- log(log(1/(1-z)))
X1 <- log(y)
Ymin <- min(Y1)-1
Ymax <- max(Y1)+0.5

#Parameter estimation
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)

#Plot
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)

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

#Maximum likelihood method

wfit1 <- fitdist(y, "weibull")
result1 <- summary(wfit1)
print(result1)
#llplot(wfit1)

#Goodness of fit
cat(" \n")
wfit2 <- fitdist(y, "weibull")
lnfit <- fitdist(y, "lnorm")
gofstat(list(wfit2, lnfit), fitnames=c("Weibull", "lognormal"))

}

 

csvデータの例

こちらをダウンロードして使ってください．