# Two step probability plot for Burr type XII distribution (BXII)
# Reference, S. Yokogawa; “Two-step probability plot for parameter estimation of lifetime distribution affected by defect clustering in time-dependent dielectric breakdown”, Japanese Journal of Applied Physics, Vol.56, pp. 07KG02-1-6 (2017).

tspbx2 <-function(s){

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

x <- read.csv("twostep.csv", header=T)
x <- c(x[,s])
y <- sort(x)

n <- length(y)
cat("Number of data")
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)
}

#2-step plot
a1 <- 1 #initial
b1 <- 1 #initial
h1 <- 1 #Initial

clm <- function(a1){
m1 <- NULL
m2 <- NULL
cu1 <- as.integer(0.000001*n)
cu2 <- as.integer(1.0*n)
m1 <- log(y)[cu1:cu2]
m2 <- log((1-z)^(-1/a1)-1)[cu1:cu2]
cor(m1,m2)
}
result1 <- optim(par=0.2,fn=clm,control=list(fnscale=-1),method="L-BFGS-B", upper=1000, lower=0.01)
a2 <- as.numeric(result1[1])

#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((1-(y0/100))^(-1/a2)-1)

#Plotting
Y1 <- log((1-z)^(-1/a2)-1)
X1 <- log(y)
xlimmin <- min(y)/10
xlimmax <- max(y)*10
ylimmin <- min(log((1-z/10)^(-1/a2)-1))
ylimmax <- max(Y1)
par(family="serif")
plot(y,Y1,log="x",xlim=c(xlimmin,xlimmax),ylim=c(ylimmin,ylimmax), xlab="Time",ylab="Cumulative probability [%]",yaxt="n")

#Grid
abline(h=probs,col="gray",lwd=1,lty=3)
axis(2,at=probs,labels=y0)

#Parameter estimation
result2 <- lsfit(Y1, X1, wt=NULL, intercept=TRUE, tolerance=1e-7, yname=NULL)
cf <- as.numeric(coef(result2))
b2 <- 1.0/cf[2]
h2 <- exp(cf[1]-1/b2*log(a2))
t1 <- h2*(a2*((1-0.001)^(-1/a2)-1))^(1/b2) #Estimation of t0.1
t2 <- h2*(a2*((1-0.000001)^(-1/a2)-1))^(1/b2) #Estimation of t1ppm

cat("Two-step estimators \n")
cat("alpha")
print(a2)
cat("beta ")
print(b2)
cat("eta ")
print(h2)
cat("t0.1 ")
print(t1)
cat("t1ppm")
print(t2)

#Line fitting
lx <- c(xlimmin,xlimmax)
Yestmin <- b2*log(xlimmin)-b2*log(h2)-log(a2)
Yestmax <- b2*log(xlimmax)-b2*log(h2)-log(a2)
ly <- c(Yestmin,Yestmax)
lines(lx,ly)

c <- cor(X1,Y1) # Correlation coefficient
cat("CC ")
print(c)
cat(" \n")

#Likelihood function of BXII
loglikelihood <- function(alpha,beta,eta){
-sum(log(beta/eta*(y/eta)^(beta-1)*(1+1/alpha*(y/eta)^beta)^(-1*alpha-1)))
}

#Maximize log-likelihood of BXII
library(bbmle)
res <- mle2(loglikelihood,start=list(alpha=a2,beta=b2,eta=h2))

rslt <- summary(res)
print(rslt)
cat(" \n")
cat("Approximate variance-covariance matrix: \n")
cv<-vcov(res)
print(cv)

}

csvデータの例

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