densities <- function (beta = c(0.5, 1, 1.5, 2, 3.6, 7), alpha = 10000) { # This function plots a variety of Weibull densities, # in particular the plot in Figure 1 of # Inference for the Weibull Distribution. # par(mar = c(0, 0, 3, 0.5) + 0.01) colors <- c("cyan", "blue", "green", "orange", "red", "purple", "black") x <- seq(0, 3 * alpha, length.out = 300) y <- dweibull(x, 1, alpha) plot(x, y, xlab = "", ylab = "", type = "l", axes = F, ylim = c(0, 3 * max(y)), main = "Weibull densities") for (i in 1:length(beta)) { y <- dweibull(x, beta[i], alpha) lines(x, y, col = colors[i]) M <- max(y) } abline(h = 0) abline(v = 0) bet <- beta legend("topright", col = colors, legend = c(expression(beta[1] == 0.5), expression(beta[2] == 1), expression(beta[3] == 1.5), expression(beta[4] == 2), expression(beta[5] == 3.6), expression(beta[6] == 7), expression(alpha == 10000)), bty = "n", lty = c(rep(1, 6), 2)) abline(v = alpha, lty = 2) text(1.01 * alpha, -5e-06, expression(alpha), adj = 0) text(1.02 * alpha, 1.1 * M, "36.8%", adj = 0) text(0.98 * alpha, 1.1 * M, "63.2%", adj = 1) } edf.plot <- function (n = 10, alpha = 10000, beta = 2) { # This function illustrates the Kolmogorov-Smirnov goodness of fit # test with samples from the Weibull distribution. Used in Figures 5 & 6 # of Inference for the Weibull Distribution. # x <- sort(rweibull(n, beta, alpha)) plot.stepfun(stepfun(x, 0:n/n), xlim = c(0, max(x, 2.5 * alpha)), pch = 16, cex = exp(-0.05 * sqrt(n)), main = "", ylab = "cumulative distribution function", cex.lab = 1.3, cex.axis = 1.3) xx <- seq(0, 3 * alpha, length.out = 1000) yy <- pweibull(xx, beta, alpha) lines(xx, yy, col = "blue", lty = 2) w.out <- Weibull.mle(x) alpha.hat <- w.out[[1]][1] beta.hat <- w.out[[1]][2] yyy <- pweibull(xx, beta.hat, alpha.hat) lines(xx, yyy, col = "orange") points(x, rep(-0.01, n), pch = 16, col = "purple", cex = 0.5) abline(h = 0) xfit <- pweibull(x, beta.hat, alpha.hat) Dm <- max(xfit - (1:n - 1)/n) Dp <- max((1:n/n - xfit)) D <- max(Dp, Dm) D.loc <- x[(1:n/n - xfit) == D | (xfit - (1:n - 1)/n) == D] y1 <- pweibull(D.loc, beta.hat, alpha.hat) if (D == Dp) { y2 <- (1:n)[x == D.loc]/n } else { y2 <- (1:n - 1)[x == D.loc]/n } segments(D.loc, y1, D.loc, y2, lwd = 2, col = "red") points(D.loc, y1, pch = 16, cex = 0.5, col = "red") points(D.loc, y2, pch = 16, cex = 0.5, col = "red") legend("topleft", legend = c(expression("EDF" == hat(F)[n]), expression("True Sampled CDF" == F[list(alpha, beta)](x)), expression("Weibull Fitted CDF" == F[list(hat(alpha), hat(beta))](x)), "KS-Distance"), col = c("black", "blue", "orange", "red"), lty = c(1, 2, 1, 1), lwd = c(1, 1, 1, 2), cex = 1.5) text(2.5 * alpha, 0.05, cex = 1.5, expression("KS-Distance" == sup({abs("" ~ hat(F)[n](x) - F[list(hat(alpha), hat(beta))](x) ~ "")}, x)), adj = 1) } GOF.AD.test <- function (AD.n = c(0.025, 0.25, 0.6, 0.9, 1.4), graphic = F) { # It is a assumed that a Weibull test of fit is # performed using the Anderson-Darling test criterion AD # and it is assumed that this criterion has been adjusted # by the factor (1+.2/sqrt(n)), i.e., AD.n = (1+.2/sqrt(n))*AD. # In this adjusted form the null distribution of AD.n is # stable in n and we have quantiles q.75, q.90, q.95, q.975, q.99 # from Stephens (1986), Tests Based on EDF Statistics, # in Goodness-of-Fit Techniques, D'Agostino and Stephens, # editors, Dekker, Table 4.17, p. 146. # For any value in AD.n the corresponding p-value is returned. # It is computed by fitting a quadratic to the quantiles against # log(p/(1-p)) where p is the upper tail probability. Using the # fitted quadratic we can read off the p correponding to any # value on the quantile axis. # The option graphic = T allows to view the interpolation # process. It was used to produce the plots in Figure 7 of # Inference for the Weibull Distribution. # Outside the range of tabulated quantiles only linear # extrapolation is done. AD.q <- c(0.474, 0.637, 0.757, 0.877, 1.038) pvec <- c(0.25, 0.1, 0.05, 0.025, 0.01) lp <- log(pvec/(1 - pvec)) coef <- lsfit(cbind(AD.q, AD.q^2), lp)$coef y0 <- rep(0, length(AD.n)) sel0.0 <- AD.n >= AD.q[1] & AD.n <= AD.q[5] sel0.1 <- AD.n < AD.q[1] sel0.2 <- AD.n > AD.q[5] y1 <- coef[1] + coef[2] * AD.q[1] + coef[3] * AD.q[1]^2 y2 <- coef[1] + coef[2] * AD.q[5] + coef[3] * AD.q[5]^2 y0[sel0.0] <- coef[1] + coef[2] * AD.n[sel0.0] + coef[3] * AD.n[sel0.0]^2 y0[sel0.1] <- y1 + (coef[2] + 2 * coef[3] * AD.q[1]) * (AD.n[sel0.1] - AD.q[1]) y0[sel0.2] <- y2 + (coef[2] + 2 * coef[3] * AD.q[5]) * (AD.n[sel0.2] - AD.q[5]) if (graphic == T) { plot(AD.q, lp, axes = F, xlab = expression(A^2 %*% (1 + 0.2/sqrt(n))), ylab = "", ylim = c(-7, 2), xlim = c(0, 1.5), pch = 16, col = "blue", cex.lab = 1.3, cex = 1.3) mtext(expression("tail probability p on " ~ log(p/(1 - p)) ~ " scale"), 2, 2.7, cex = 1.3) axis(1, cex.axis = 1.3) pvec.lab <- c(0.9, 0.75, 0.5, pvec, 0.001) lp.lab <- log(pvec.lab/(1 - pvec.lab)) axis(2, at = lp.lab, labels = pvec.lab, cex.axis = 1.3) x <- seq(0, 1.5, 0.001) y <- coef[1] + coef[2] * x + coef[3] * x^2 sel0 <- x >= AD.q[1] & x <= AD.q[5] sel1 <- x < AD.q[1] sel2 <- x > AD.q[5] ye1 <- y1 + (coef[2] + 2 * coef[3] * AD.q[1]) * (x - AD.q[1]) ye2 <- y2 + (coef[2] + 2 * coef[3] * AD.q[5]) * (x - AD.q[5]) lines(x[sel1], ye1[sel1], lty = 2, col = "blue") lines(x[sel2], ye2[sel2], lty = 2, col = "blue") lines(x[sel0], y[sel0], col = "blue") lines(x[sel1], y[sel1], col = "blue", lty = 3) lines(x[sel2], y[sel2], col = "blue", lty = 3) points(AD.n, y0, col = "red", pch = 16, cex = 1.3) legend("topright", legend = c("tabled values", "interpolated/extrapolated values"), pch = c(16, 16), col = c("blue", "red"), bty = "n", cex = 1.4) } pvalues <- exp(y0)/(1 + exp(y0)) pvalues } GOF.CvM.test <- function (CvM.n = c(0.025, 0.07, 0.1, 0.15, 0.25), graphic = F) { # It is a assumed that a Weibull test of fit is # performed using the Cramer-von Mises test criterion CvM # and it is assumed that this criterion has been adjusted # by the factor (1+.2/sqrt(n)), i.e., CvM.n = (1+.2/sqrt(n))*CvM. # In this adjusted form the null distribution of CvM.n is # stable in n and we have quantiles q.75, q.90, q.95, q.975, q.99 # from Stephens (1986), Tests Based on EDF Statistics, # in Goodness-of-Fit Techniques, D'Agostino and Stephens, # editors, Dekker, Table 4.17, p. 146. # For any value in CcM.n the corresponding p-value is returned. # It is computed by fitting a quadratic to the quantiles against # log(p/(1-p)) where p is the upper tail probability. Using the # fitted quadratic we can read off the p correponding to any # value on the quantile axis. # The option graphic = T allows to view the interpolation # process. It was used to produce the plots in Figure 7 # of Inference for the Weibull Distribution. # Outside the range of tabulated quantiles only linear # extrapolation is done. CvM.q <- c(0.073, 0.102, 0.124, 0.146, 0.175) pvec <- c(0.25, 0.1, 0.05, 0.025, 0.01) lp <- log(pvec/(1 - pvec)) coef <- lsfit(cbind(CvM.q, CvM.q^2), lp)$coef y1 <- coef[1] + coef[2] * CvM.q[1] + coef[3] * CvM.q[1]^2 y2 <- coef[1] + coef[2] * CvM.q[5] + coef[3] * CvM.q[5]^2 sel0.0 <- CvM.n >= CvM.q[1] & CvM.n <= CvM.q[5] sel0.1 <- CvM.n < CvM.q[1] sel0.2 <- CvM.n > CvM.q[5] y0 <- rep(0, length(CvM.n)) y0[sel0.0] <- coef[1] + coef[2] * CvM.n[sel0.0] + coef[3] * CvM.n[sel0.0]^2 y0[sel0.1] <- y1 + (coef[2] + 2 * coef[3] * CvM.q[1]) * (CvM.n[sel0.1] - CvM.q[1]) y0[sel0.2] <- y2 + (coef[2] + 2 * coef[3] * CvM.q[5]) * (CvM.n[sel0.2] - CvM.q[5]) if (graphic == T) { plot(CvM.q, lp, axes = F, xlab = expression(W^2 %*% (1 + 0.2/sqrt(n))), ylab = "", ylim = c(-7, 2), xlim = c(0, 0.3), pch = 16, col = "blue", cex.lab = 1.3, cex = 1.3) mtext(expression("tail probability p on " ~ log(p/(1 - p)) ~ " scale"), 2, 2.7, cex = 1.3) axis(1, cex.axis = 1.3) pvec.lab <- c(0.9, 0.75, 0.5, pvec, 0.001) lp.lab <- log(pvec.lab/(1 - pvec.lab)) axis(2, at = lp.lab, labels = pvec.lab, cex.axis = 1.3) x <- seq(0, 0.3, 0.001) y <- coef[1] + coef[2] * x + coef[3] * x^2 sel0 <- x >= CvM.q[1] & x <= CvM.q[5] sel1 <- x < CvM.q[1] sel2 <- x > CvM.q[5] ye1 <- y1 + (coef[2] + 2 * coef[3] * CvM.q[1]) * (x - CvM.q[1]) ye2 <- y2 + (coef[2] + 2 * coef[3] * CvM.q[5]) * (x - CvM.q[5]) lines(x[sel1], ye1[sel1], lty = 2, col = "blue") lines(x[sel2], ye2[sel2], lty = 2, col = "blue") lines(x[sel0], y[sel0], col = "blue") lines(x[sel1], y[sel1], col = "blue", lty = 3) lines(x[sel2], y[sel2], col = "blue", lty = 3) points(CvM.n, y0, col = "red", pch = 16, cex = 1.3) legend("topright", legend = c("tabled values", "interpolated/extrapolated values"), pch = c(16, 16), col = c("blue", "red"), bty = "n", cex = 1.3) } pvalues <- exp(y0)/(1 + exp(y0)) pvalues } GOF.KS.test <- function (D = c(0.7, 0.75, 0.85, 0.95, 1, 1.05), n = 25, graphic = F) { # It is assumed the a Weibull test of fit is performed # using the Kolmogorov-Smirnov test of fit metric KS. # This function calculates p-values for the sqrt(n) adjusted # value of KS, i.e., it is assumed that D = sqrt(n)*KS. # In spite of the sqrt(n) factor the null distribution # of this statistic still varies with n. # A two-fold interpolation and/or extrapolation is needed in # order to calculate the p-values for the entered values of D. # First one interpolates the quantiles for the given sample size # n, doing a quadratic interpolation of the given quantiles # versus 1/sqrt(n). Then one does a cubic interpolation of the # log-odds of the upper tail probability against these # interpolated quantiles in order to arrive at the desired # p-values. Outside the range of the quantiles only linear # extraplation is used. # The tabled quantiles q.90, q.95, q.975, and q.99 for # n=10,20,50,Inf are taken from the Stephens (1986), # Tests Based on EDF Statistics, in Goodness-of-Fit Techniques, # D'Agostino and Stephens editors, Dekker, Table 4.18, p. 148. # The returned value is the vector of p-values corresponding to # the sqrt(n) adjusted KS vector D and the given sample size # n. # The option graphic = T allows to see the interpolation # plots. It was used to produce the plots in Figure 8 of # Inference for the Weibull Distribution. # D.tab <- rbind(c(0.76, 0.819, 0.88, 0.944), c(0.779, 0.843, 0.907, 0.973), c(0.79, 0.856, 0.922, 0.988), c(0.803, 0.874, 0.939, 1.007)) nvec <- c(10, 20, 50, Inf) pvec <- c(0.1, 0.05, 0.025, 0.01) lp <- log(pvec/(1 - pvec)) lp0 <- rep(0, length(D)) xn <- 1/sqrt(nvec) xn2 <- 1/nvec xn0 <- 1/sqrt(n) xn02 <- 1/n coef.mat <- NULL coefp.mat <- NULL q0 <- NULL for (i in 1:4) { coef <- lsfit(cbind(xn, xn2), D.tab[, i])$coef coef.mat <- rbind(coef.mat, coef) if (n >= 10) { q0 <- c(q0, coef[1] + coef[2] * xn0 + coef[3] * xn02) } else { x01 <- 1/sqrt(10) y01 <- coef[1] + coef[2] * x01 + coef[3] * x01^2 q0 <- c(q0, y01 + (coef[2] + 2 * coef[3] * x01) * (xn0 - x01)) } } coefp <- lsfit(cbind(q0, q0^2, q0^3), lp)$coef lp01 <- coefp[1] + coefp[2] * q0[1] + coefp[3] * q0[1]^2 + coefp[4] * q0[1]^3 lp02 <- coefp[1] + coefp[2] * q0[4] + coefp[3] * q0[4]^2 + coefp[4] * q0[4]^3 sel.int <- D >= q0[1] & D <= q0[4] sel.ext1 <- D < q0[1] sel.ext2 <- D > q0[4] lp0[sel.int] <- coefp[1] + coefp[2] * D[sel.int] + coefp[3] * D[sel.int]^2 + coefp[4] * D[sel.int]^3 lp0[sel.ext1] <- lp01 + (coefp[2] + 2 * coefp[3] * q0[1] + 3 * coefp[4] * q0[1]^2) * (D[sel.ext1] - q0[1]) lp0[sel.ext2] <- lp02 + (coefp[2] + 2 * coefp[3] * q0[4] + 3 * coefp[4] * q0[4]^2) * (D[sel.ext2] - q0[4]) pvalues <- exp(lp0)/(1 + exp(lp0)) if (graphic == T) { for (i in 1:4) { Di <- D.tab[i, ] Di2 <- Di^2 Di3 <- Di^3 coefp.mat <- rbind(coefp.mat, lsfit(cbind(Di, Di2, Di3), lp)$coef) } plot(xn, D.tab[, 1], ylim = c(min(0.6, D.tab), max(D.tab)), xlim = c(0, 0.5), ylab = "", xlab = expression(n ~ " (on " ~ 1/sqrt(n) ~ " scale)"), axes = F, main = expression("quadratic interpolation & linear extrapolation in " ~ 1/sqrt(n)), pch = 16, col = "blue", cex.main = 1.5, cex.lab = 1.3, cex = 1.3) axis(1, pos = min(0.6, D.tab), at = 1/sqrt(c(4:10, 15, 20, 25, 30, 40, 50, 100, 200, 500, Inf)), labels = c(4:10, 15, 20, 25, 30, 40, 50, 100, 200, 500, expression(infinity)), cex.axis = 1.3) axis(2, pos = -0.005, at = D.tab[4, ], labels = D.tab[4, ], cex.axis = 1.3) mtext(expression(sqrt(n) %*% ~"D-quantile"), 2, 2.5, cex = 1.5) mtext(c(expression(D[0.9]), expression(D[0.95]), expression(D[0.975]), expression(D[0.99])), 2, 1.2, at = D.tab[4, ], cex = 1.3) abline(v = -0.005) abline(h = min(0.6, D.tab)) xx <- seq(0, 0.5, 0.001) for (i in 1:4) { yy <- coef.mat[i, 1] + coef.mat[i, 2] * xx + coef.mat[i, 3] * xx^2 lines(xx[xx <= 1/sqrt(10)], yy[xx <= 1/sqrt(10)], col = "blue") yy2 <- coef.mat[i, 1] + coef.mat[i, 2]/sqrt(10) + coef.mat[i, 3]/10 yye <- yy2 + (coef.mat[i, 2] + 2 * coef.mat[i, 3]/sqrt(10)) * (xx[xx > 1/sqrt(10)] - 1/sqrt(10)) lines(xx[xx > 1/sqrt(10)], yye, lty = 2, col = "blue") lines(xx[xx > 1/sqrt(10)], yy[xx > 1/sqrt(10)], col = "blue", lty = 3) points(xn, D.tab[, i], pch = 16, col = "blue", cex = 1.3) } legend("topright", legend = c("tabled values", "interpolated/extrapolated values"), pch = c(16, 16), col = c("blue", "red"), bty = "n", cex = 1.3) points(rep(xn0, 4), q0, pch = 16, col = "red", cex = 1.3) readline("hit return\n") plot(D.tab[1, ], lp, xlab = expression(sqrt(n) %*% D), ylab = "", xlim = c(0.7, 1.05), ylim = c(-7, -1), axes = F, main = "cubic interpolation & linear extrapolation in D", cex = 1.3, cex.lab = 1.3, cex.main = 1.5) axis(1, cex.axis = 1.3) mtext(expression("tail probability p on " ~ log(p/(1 - p)) ~ " scale"), 2, 2.5, cex = 1.3) pos <- c(0.2, 0.1, 0.05, 0.025, 0.01, 0.005, 0.001) axis(2, at = log(pos/(1 - pos)), labels = c(0.2, 0.1, 0.05, 0.025, 0.01, 0.005, 0.001), pos = 0.7, cex.axis = 1.3) abline(v = 0.7) for (i in 1:4) { points(D.tab[i, ], lp, pch = 16, col = "blue", cex = 1.3) } xx0 <- seq(0.7, 1.05, 1e-04) points(q0, lp, col = "red", pch = 16, cex = 1.3) for (i in 1:4) { lpxx0 <- coefp.mat[i, 1] + coefp.mat[i, 2] * xx0 + coefp.mat[i, 3] * xx0^2 + coefp.mat[i, 4] * xx0^3 lines(xx0[lpxx0 >= lp[4] & lpxx0 <= lp[1]], lpxx0[lpxx0 >= lp[4] & lpxx0 <= lp[1]], col = "blue") lp2 <- coefp.mat[i, 1] + coefp.mat[i, 2] * D.tab[i, 4] + coefp.mat[i, 3] * D.tab[i, 4]^2 + coefp.mat[i, 4] * D.tab[i, 4]^3 lp1 <- coefp.mat[i, 1] + coefp.mat[i, 2] * D.tab[i, 1] + coefp.mat[i, 3] * D.tab[i, 1]^2 + coefp.mat[i, 4] * D.tab[i, 1]^3 lpe1 <- lp1 + (coefp.mat[i, 2] + 2 * D.tab[i, 1] * coefp.mat[i, 3] + 3 * D.tab[i, 1]^2 * coefp.mat[i, 4]) * (xx0[lpxx0 > lp[1]] - D.tab[i, 1]) lpe2 <- lp2 + (coefp.mat[i, 2] + 2 * D.tab[i, 4] * coefp.mat[i, 3] + 3 * D.tab[i, 4]^2 * coefp.mat[i, 4]) * (xx0[lpxx0 < lp[4]] - D.tab[i, 4]) lines(xx0[lpxx0 > lp[1]], lpe1, col = "blue", lty = 2) lines(xx0[lpxx0 < lp[4]], lpe2, col = "blue", lty = 2) } lpq0 <- coefp[1] + coefp[2] * xx0 + coefp[3] * xx0^2 + coefp[4] * xx0^3 lines(xx0[lpq0 >= lp[4] & lpq0 <= lp[1]], lpq0[lpq0 >= lp[4] & lpq0 <= lp[1]], col = "red") lpn1 <- coefp[1] + coefp[2] * q0[1] + coefp[3] * q0[1]^2 + coefp[4] * q0[1]^3 lpn2 <- coefp[1] + coefp[2] * q0[4] + coefp[3] * q0[4]^2 + coefp[4] * q0[4]^3 lpne1 <- lpn1 + (coefp[2] + 2 * q0[1] * coefp[3] + 3 * q0[1]^2 * coefp[4]) * (xx0[xx0 < q0[1]] - q0[1]) lines(xx0[xx0 < q0[1]], lpne1, col = "red", lty = 2) lpne2 <- lpn2 + (coefp[2] + 2 * q0[4] * coefp[3] + 3 * q0[4]^2 * coefp[4]) * (xx0[xx0 > q0[4]] - q0[4]) lines(xx0[xx0 > q0[4]], lpne2, col = "red", lty = 2) legend("topright", legend = c("tabled values", "interpolated quantiles", "interpolated/extrapolated values"), pch = c(16, 16, 16), col = c("blue", "red", "black"), bty = "n", cex = 1.3) points(D, lp0, pch = 16, col = "black", cex = 1.3) } pvalues } kazbek <- structure(list(time = c(186000L, 166000L, 159000L, 255000L, 255000L, 255000L, 255000L, 255000L, 255000L, 255000L, 255000L, 255000L, 255000L), status = c(1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("time", "status"), class = "data.frame", row.names = c(NA, -13L)) # this data set gives right censored data for which the classical # MLE confidence bounds are not monotone. KS.CvM.AD <- function (x) { # this function computes the three test of fit criteria # KS, CvM and AD # if (!exists("survreg")) library(survival) xs <- sort(x) n <- length(xs) statusx <- rep(1, n) dat.weibull <- data.frame(xs, statusx) names(dat.weibull) <- c("time", "status") out.weibull <- survreg(Surv(time, status) ~ 1, dist = "weibull", data = dat.weibull) u.hat <- out.weibull$coef b.hat <- out.weibull$scale V <- (log(xs) - u.hat)/b.hat Z <- 1 - exp(-exp(V)) AD <- -n - sum((2 * (1:n) - 1) * (log(Z) + log(1 - rev(Z))))/n CvM <- sum((Z - (2 * (1:n) - 1)/(2 * n))^2) + 1/(12 * n) KS <- max(max((1:n)/n - Z), max(Z - ((1:n) - 1)/n)) c(KS, CvM, AD) } shockabsorber <- structure(list(time = c(6700L, 6950L, 7820L, 8790L, 9120L, 9660L, 9820L, 11310L, 11690L, 11850L, 11880L, 12140L, 12200L, 12870L, 13150L, 13330L, 13470L, 14040L, 14300L, 17520L, 17540L, 17890L, 18450L, 18960L, 18980L, 19410L, 20100L, 20100L, 20150L, 20320L, 20900L, 22700L, 23490L, 26510L, 27410L, 27490L, 27890L, 28100L ), status = c(1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L)), .Names = c("time", "status"), row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38"), class = "data.frame") # This data set represents right censored data and is taken # from Meeker and Escobar, Statistical Methods for Reliability Data, # John Wiley & Sons, New York. Table 3.4, p. 59 test40 <- function () { # This function produces the plot in Figure 9 of # Inference for the Weibull Distribution. # par(mfrow = c(2, 1)) x <- c(1.628, 1.886, 2.246, 2.486, 2.796, 3.2, 3.761, 4.565, 4.826) p0 <- 0.96 p1 <- 0.98 x0 <- x x0[9] <- 5.725 p <- c(0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.98) plot(x, log(p/(1 - p)), xlim = c(min(x, x0), max(x, x0)), xlab = expression("Bain's tabled quantiles for n=40," ~ gamma == 0.9), cex = 1.35, pch = 1, col = "red") points(4.826, log(p0/(1 - p0)), pch = 16, col = "red", cex = 1.35) text(4.826 + 0.05, log(p0/(1 - p0)), p0, adj = 0) text(4.826 + 0.05, log(p1/(1 - p1)), p1, adj = 0) plot(x0, log(p/(1 - p)), xlim = c(min(x, x0), max(x, x0)), xlab = expression("our simulated quantiles for n=40," ~ gamma == 0.9), pch = 16, col = "blue") points(x, log(p/(1 - p)), pch = 1, cex = 1.35, col = "red") } timing.plot <- function () { # This function produces the plot in Figure 4 of # Inference for the Weibull Distribution. # x <- c(10, 100, 500, 1000, 5000) y <- c(1.53, 1.71, 3.47, 7.15, 26.68)/1000 plot(x, y, xlab = "sample size n", ylab = "time to compute Weibull mle's (sec)", pch = 16, col = "purple", ylim = c(0, 0.03), cex.lab = 1.5, cex.axis = 1.5, cex = 1.5) ls.out <- lsfit(x, y) abline(ls.out, col = "blue") text(0, 0.029, paste("intercept = ", format(signif(ls.out$coef[1], 4)), ", \n\t\tslope = ", format(signif(ls.out$coef[2], 4))), adj = 0, cex = 1.5) ls.out$coef } WEIB.base <- function (dat = weibproblem.dat, gam = 0.95, p = 0.1, q0 = NULL, q1 = NULL, q2 = NULL, qr = NULL, comp.mode = 1, input = NULL) { # This routine uses survreg from the package survival. # Also needed to reference the shape parameter differently. # gam = .95 gives the confidence level for the intervals # p gives the probability level for the quantile of interest # q0 gives the threshold for left tail probability calculation # q1 < q2 give the thresholds for the calculation of the # probability that life will end by time q2 when it is known to # have exceeded q1. # qr gives the life times at which failure rates are # calculated", # Note that the following relation holds between failure rates # of Weibull and extreme value models: # lambda_W(t) = lambda_E[log(t)]/t # comp.mode = 1 indicates that survreg is to be run to # calculate the m.l.e.'s and accompanying items. # comp.mode = 2 indicates that survreg should not be run and # that the information from survreg is supplied via the input # argument, giving it the $out component from the output of a # previous run of WEIB.base as input. # This function is used only in weibull.paper.mon to construct # bounds for p-quantiles or tail probabilities by three different # methods if (!exists("survreg")) library(survival) interval.xp <- NULL interval.p <- NULL xp <- NULL phat <- NULL cond.prob <- NULL interval.cond.prob <- NULL lambda <- NULL interval.lambda <- NULL names(gam) <- "confidence level" attach(dat) n <- length(time) wp <- log(-log(1 - p)) z <- qnorm((1 - gam)/2) zvec <- c(-z, z) if (comp.mode == 1) { out <- survreg(Surv(time, status) ~ 1, dist = "weibull", data = dat) } else { out <- input } beta <- 1/out$scale alpha <- exp(out$coef) names(alpha) <- c("alpha (characteristic life)") names(beta) <- c("beta (shape)") interval.b <- 1/exp(log(out$scale) + zvec * sqrt(out$var[2, 2])) interval.a <- exp(out$coef - zvec * sqrt(out$var[1, 1])) names(interval.a) <- c("alpha.L", "alpha.U") names(interval.b) <- c("beta.L", "beta.U") mu <- out$coef b <- out$scale if (length(p) > 0) { zp <- rep(z, length(p)) zvecp <- cbind(zp, -zp) xp <- mu + b * wp interval.xp <- xp + zvecp * sqrt(out$var[1, 1] + 2 * out$var[1, 2] * wp * b + out$var[2, 2] * wp^2 * b^2) interval.xp <- exp(interval.xp) interval.xp <- cbind(p, interval.xp) dimnames(interval.xp) <- list(NULL, c("p", "xp.L", "xp.U")) xp <- cbind(p, exp(xp)) dimnames(xp) <- list(NULL, c("p", "xp")) } if (length(q0) > 0) { zp <- rep(z, length(q0)) zvecp <- cbind(zp, -zp) x0 <- log(q0) phat <- (x0 - mu)/b interval.p <- phat + zvecp * (1/b) * sqrt(out$var[1, 1] + 2 * out$var[1, 2] * (x0 - mu) + (x0 - mu)^2 * out$var[2, 2]) phat <- 1 - exp(-exp(phat)) interval.p <- 1 - exp(-exp(interval.p)) interval.p <- cbind(q0, interval.p) dimnames(interval.p) <- list(NULL, c("q0", "p.L", "p.U")) phat <- cbind(q0, phat) dimnames(phat) <- list(NULL, c("q0", "prob")) } if (length(q1) == length(q2) & length(q1) > 0) { zp <- rep(z, length(q1)) zvecp <- cbind(zp, -zp) theta <- (log(q1) - mu)/b + log(exp(log(q2/q1)/b) - 1) cond.prob <- 1 - exp(-exp(theta)) cond.prob <- cbind(q1, q2, cond.prob) dimnames(cond.prob) <- list(NULL, c("q1", "q2", "cond.prob")) interval.cond.prob <- theta + zvecp * (1/b) * sqrt(out$var[1, 1] + 2 * out$var[1, 2] * (log(q1) - mu + (exp(log(q2/q1)/b) * log(q2/q1))/(exp(log(q2/q1)/b) - 1)) + out$var[2, 2] * (log(q1) - mu + (exp(log(q2/q1)/b) * log(q2/q1))/(exp(log(q2/q1)/b) - 1))^2) interval.cond.prob <- 1 - exp(-exp(interval.cond.prob)) interval.cond.prob <- cbind(q1, q2, interval.cond.prob) dimnames(interval.cond.prob) <- list(NULL, c("q1", "q2", "cond.prob.L", "cond.prob.U")) } if (length(qr) > 0) { zp <- rep(z, length(qr)) zvecp <- cbind(zp, -zp) f.rate <- (-log(b) + (1/b - 1) * log(qr) - mu/b) interval.f.rate <- f.rate + zvecp * sqrt(out$var[1, 1]/b^2 + 2 * out$var[1, 2] * (1/b) * (1 + (1/b) * log(qr) - mu/b) + out$var[2, 2] * (1 + (1/b) * log(qr) - mu/b)^2) lambda <- cbind(qr, exp(f.rate)) interval.lambda <- cbind(qr, exp(interval.f.rate)) dimnames(lambda) <- list(NULL, c("qr", "lambda")) dimnames(interval.lambda) <- list(NULL, c("qr", "lambda.L", "lambda.U")) } result <- list(gam = gam, alpha = alpha, beta = beta, interval.alpha = interval.a, interval.beta = interval.b, xp = xp, interval.xp = interval.xp, prob = phat, interval.p = interval.p, cond.prob = cond.prob, interval.cond.prob = interval.cond.prob, lambda = lambda, interval.lambda = interval.lambda, out = out) detach() result } weib.sample <- c(13557.9878403624, 13186.3684747875, 11991.4191645563, 8145.2536210027, 2116.57665045492, 16113.6417181945, 9927.16122041184, 11221.9776364315, 9658.0979793509, 4535.92555829378) # uncensored Weibull data set weibproblem.dat <- structure(list(time = c(626.1, 651.7, 684.7, 686.3, 698.2, 707.7, 709.8, 714.7, 718, 719.6, 720.9, 721.9, 726.7, 740.3, 752.9, 760.3, 764, 764.8, 768.3, 773.6, 774.4, 784.1, 785.3, 788.9, 790.3, 793.2, 794, 806.1, 816.2, 825, 826.5, 829.8, 832.3, 839.4, 840.5, 843.1, 845.2, 849.1, 849.2, 856.2, 856.8, 859.1, 868.9, 869.3, 881.1, 887.8, 890.5, 898.2, 921.5, 934.8), status = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("time", "status"), class = "data.frame", row.names = c(NA, -50L)) # right censored Weibull data set Weibull.GOF.test <- function (x) { # This test perform a Weibull goodness of fit test # for the provided sample x. It does this using the # Kolomogorov-Smirnov, Cramer-von Mises and Anderson-Darling # test of fit criteria, which are adjusted appropriately # for n prior to their p-value calculation via GOF.KS.test, # GOF.CvM.test, and GOF.AD.test, which are assumed to be in the # work space. if (!exists("survreg")) library(survival) xs <- sort(x) n <- length(xs) statusx <- rep(1, n) dat.weibull <- data.frame(xs, statusx) names(dat.weibull) <- c("time", "status") out.weibull <- survreg(Surv(time, status) ~ 1, dist = "weibull", data = dat.weibull) u.hat <- out.weibull$coef b.hat <- out.weibull$scale xs <- rev(xs) Z <- exp(-exp(-(-log(xs) + u.hat)/b.hat)) AD <- -n - sum((2 * (1:n) - 1) * (log(Z) + log(1 - rev(Z))))/n CvM <- sum((Z - (2 * (1:n) - 1)/(2 * n))^2) + 1/(12 * n) KS <- max(max((1:n)/n - Z), max(Z - ((1:n) - 1)/n)) p.AD <- GOF.AD.test(AD * (1 + 0.2/sqrt(n))) p.CvM <- GOF.CvM.test(CvM * (1 + 0.2/sqrt(n))) p.KS <- GOF.KS.test(sqrt(n) * KS, n) pvalues <- c(p.KS, p.CvM, p.AD) names(pvalues) <- c("pval.KS", "pval.CvM", "pval.AD") pvalues } Weibull.GOF.test.sim <- function (n = 10, Nsim = 1000) { # This function simulates samples of size n # from a Weibull distribution with shape 1 # and scale 1 and for each such sample it computes the # p-value as produced by Weibull.GOF.test. # It does this Nsim times and then shows the respective # histograms of these p-values. At least at the low end # (say [0,.25]) one expects to see roughly a uniform # distribution. pval <- matrix(rep(0, Nsim * 3), ncol = 3) for (i in 1:Nsim) { pval[i, ] <- Weibull.GOF.test(rweibull(n, 1)) } par(mfrow = c(3, 1)) hist(pval[, 1], nclass = 50, xlim = c(0, 1), xlab = paste("p-values of KS test ( n =", n, ")"), main = paste("Weibull KS Test of Fit, Nsim =", Nsim, " Simulations"), col = c("blue", "orange")) hist(pval[, 2], nclass = 50, xlim = c(0, 1), xlab = paste("p-values of CvM test ( n =", n, ")"), main = paste("Weibull CvM Test of Fit, Nsim =", Nsim, " Simulations"), col = c("blue", "orange")) hist(pval[, 3], nclass = 50, xlim = c(0, 1), xlab = paste("p-values of AD test ( n =", n, ")"), main = paste("Weibull AD Test of Fit, Nsim =", Nsim, " Simulations"), col = c("blue", "orange")) } Weibull.mle <- function (x = NULL, n = NULL) { # This function computes the maximum likelihood estimates of # alpha and beta for complete or type II censored samples # assumed to come from a 2-parameter Weibull distribution. # Here x is the sample, either the full sample or the first # r observations of a type II censored sample. In the latter # case one must specify the full sample size n, otherwise x is # treated as a full sample. # If x is not given then a default full sample of size n = 10, # namely # c(7,12.1,22.8,23.1,25.7,26.7,29.0,29.9,39.5,41.9) is analyzed # and the returned # results should be # $mles # alpha.hat beta.hat # 28.914017 2.799793 # In the type II censored usage # Weibull.mle(c(7,12.1,22.8,23.1,25.7),10) # $mles # alpha.hat beta.hat # 30.725992 2.432647 # Fritz Scholz, May 2008. if (is.null(x)) x <- c(7, 12.1, 22.8, 23.1, 25.7, 26.7, 29, 29.9, 39.5, 41.9) r <- length(x) if (is.null(n)) { n <- r } else { if (r > n || r < 2) { return("x must have length r with: 2 <= r <= n") } } xs <- sort(x) if (!exists("survreg")) library(survival) if (r < n) { statusx <- c(rep(1, r), rep(0, n - r)) dat.weibull <- data.frame(c(xs, rep(xs[r], n - r)), statusx) } else { statusx <- rep(1, n) dat.weibull <- data.frame(xs, statusx) } names(dat.weibull) <- c("time", "status") out.weibull <- survreg(Surv(time, status) ~ 1, dist = "weibull", data = dat.weibull) alpha.hat <- exp(out.weibull$coef) beta.hat <- 1/out.weibull$scale parms <- c(alpha.hat, beta.hat) names(parms) <- c("alpha.hat", "beta.hat") list(mles = parms) } Weibull.monotone.R <- function (dat, gam = 0.95, cov0 = NULL, p = 0.1, x = NULL, intercept = T) { # This routine does a Weibull maximum likelihood analysis on # right censored positive lifetime data with or without # additional covariates. It computes estimates and confidence # intervals for the Weibull shape parameter, for the regression # coefficients coef of log(alpha) = t(cov)%*%coef in relation to # the supplied covariates cov, for the mean, specified # p-quantiles and left tail probabilities # P( X < x ). # For the latter three quantities one must specify a taget # covariate vector cov0 in case additional covariates are # used. The quantile and tail probability intervals done by # classical methods are not necessarily monotone in p or x and # thus quantile and tail probability interval end points are not # inverse to each other. Intervals are also given by an # alternate method that guarantees monotone interval bounds so # that quantile lower(upper) bounds invert to tail probability # upper(lower) bounds. # This routine excecutes library(survival) if it has not yet # been done prior to the call. # INPUT: dat = # the data frame with first column named "time" # containing positive lifetimes, and with last column # named "status, giving the censoring status indicator # (0 for right censored and 1 for failure). # The columns between time and status represent # additional covariates, if present. # gam = .95 (default) # gives the confidence level for the intervals. # cov0 = NULL is a covariate vector. Its length should # match the number of covariates provided in the # data set. It will be preceded by 1 internally # if intercept = T. # p = .1 (default) gives the probability level for # the p-quantile of interest (p can be a # vector). # x = NULL # gives the threshold for the left tail probability # of interest [ P(X < x) ](x can be a vector) # # OUTPUT: stats = a vector with components # "confidence level", "n observations", # "r failures", # coef = a vector with components coef.hat # (estimate of coef) and its confidence # bounds coef.L and coef.U. # These coefficients explain log(alpha) as linear # combination of the covariates. # beta = a vector with components beta.hat # (estimate of Weibull shape parameter beta), # beta.L and beta.U (confidence interval for # beta). # mu = a vector with components mu.hat (estimate # of the Weibull mean = alpha*gamma(1+ # 1/beta) ), mu.L and mu.U (confidence # interval for mu).", # xp = a matrix with first column p, and columns # 2,...,8 containing the corresponding # quantile estimates xp.hat, classical # confidence intervals xp.L & xp.U, and the # monotone confidence intervals xp.L.m & # xp.U.m. # tail.prob = a matrix with first column x, and columns # 2,...,6 containing the corresponding # estimates of the left tail probability, # the classical confidence intervals # p.L & p.U, and the monotone confidence # intervals p.L.m & p.U.m. # time = processing time ", # The methodolgy for the classical intervals is based on the # approximate normality of coef.hat and log(beta.hat) as far as # bounds for coef[i] beta and xp are concerned. The classical # bounds for P(X < x) are either based on the approximate # normality [log(x) - log(alpha.hat)]*beta.hat or on the # approximate normality of the logit of the estimate for # P(X < x). # The methodology for the monotone bounds is explained in # Unified Confidence Bounds for Censored Weibull # Data With Covariates, # see http://faculty.washington.edu/fscholz/Reports/confbd.pdf # The reason for calling this routine: Weibull.unified is # the Weibull quantile xp = alpha(coef)* (-log(1-p))^(1/beta) or # log(xp) = u(coef)+ b*wp with wp = log(-log(1-p)), b=1/beta # and u(coef)= log(alpha(coef)) = t(cov0)%*%coef the location # parameter of the log(Weibull) data at the covariate vector # cov0. # This log(xp) can become equal to b or coef[i] by choosing # cov0=(0,...,0) and wp=1, i.e., p=1-exp(-exp(1)) # or by choosing cov0=(0,..,0,1,0,..,0) with 1 in the i-th # position and then taking wp=0, i.e., p=1-exp(-1). # This way the classical confidence bounds for beta and # coef[i] coincide with the ones obtained via the special case # xp bounds. # Fritz Scholz, July 2005 (revised August 2015). if (!exists("survreg")) library(survival) now <- proc.time()[1:2] n.time <- dim(dat)[1] ncol <- ncol(dat) r <- sum(dat[, ncol]) n.cov <- ncol - 2 names(gam) <- "confidence level" wp <- log(-log(1 - p)) interval.xp <- NULL interval.p <- NULL xp <- NULL phat <- NULL tail.prob <- NULL if (n.cov != length(cov0)) { return(paste("length of cov0 =", length(cov0), "should match", n.cov, ", the number of covariates provided in the data matrix")) } one <- 1 names(one) <- "intercept" if (n.cov > 0 & intercept == T) { names(cov0) <- names(dat)[2:(ncol - 1)] cov0 <- c(one, cov0) } if (n.cov == 0 & intercept == T) { cov0 <- one } if (n.cov > 0 & intercept == F) { names(cov0) <- names(dat)[2:(ncol - 1)] } if (n.cov > 0) { form1 <- as.formula(paste("Surv(", names(dat)[1], ",", names(dat)[ncol], ") ~", paste(names(dat)[2:(ncol - 1)], collapse = "+"), collapse = "")) form2 <- as.formula(paste("Surv(", names(dat)[1], ",", names(dat)[ncol], ") ~", paste(names(dat)[2:(ncol - 1)], collapse = "+"), "-1", collapse = "")) } if (n.cov > 0 & intercept == T) { cov <- dat[, 2:(n.cov + 1)] out <- lsfit(cov, rep(1, n.time), int = F) if (sum(out$resid^2) < 1e-07) { return("matrix of given covariates is collinear with intercept") } } if (intercept == F & length(cov0) != n.cov) { return(paste("covariate vector", paste("(", paste(cov0, collapse = ","), ")", collapse = ""), "not of correct length", n.cov)) } if (intercept == T) { if (n.cov > 0) { if (exists("is.R")) { out <- survreg(form1, dist = "weibull", dat = dat) } else { out <- survReg(form1, dist = "weibull", data = dat) } } else { if (exists("is.R")) { out <- survreg(Surv(time, status) ~ 1, dist = "weibull", data = dat) } else { out <- survReg(Surv(time, status) ~ 1, dist = "weibull", data = dat) } } } else { if (n.cov > 0) { if (exists("is.R")) { out <- survreg(form2, dist = "weibull", dat = dat) } else { out <- survReg(form2, dist = "weibull", data = dat) } } else { return("can't have intercept = F and no covariates") } } if (intercept == T) { n.cov <- n.cov + 1 } z <- qnorm((1 - gam)/2) zvec <- c(-z, z) tau.11 <- as.matrix(out$var[1:n.cov, 1:n.cov]) tau.12 <- as.vector(out$var[1:n.cov, n.cov + 1]) tau.21 <- t(as.vector(out$var[1:n.cov, n.cov + 1])) tau.22 <- out$var[n.cov + 1, n.cov + 1] b.hat <- out$scale coef.hat <- out$coef u.hat <- t(cov0) %*% coef.hat beta.hat <- 1/b.hat interval.beta <- 1/exp(log(b.hat) + zvec * sqrt(tau.22)) interval.coef <- cbind(coef.hat + z * sqrt(diag(tau.11)), coef.hat - z * sqrt(diag(tau.11))) coef <- signif(cbind(coef.hat, interval.coef), 5) beta <- signif(c(beta.hat, interval.beta), 5) dimnames(coef) <- list(NULL, c("coef.hat", "coef.L", "coef.U")) names(beta) <- c("beta.hat", "beta.L", "beta.U") theta.hat <- u.hat + log(gamma(1 + b.hat)) mu.hat <- exp(theta.hat) sig.theta.hat <- sqrt(t(cov0) %*% tau.11 %*% cov0 + 2 * b.hat * digamma(1 + b.hat) * t(cov0) %*% tau.12 + b.hat^2 * (digamma(1 + b.hat))^2 * tau.22) theta.int <- theta.hat - zvec * sig.theta.hat mu <- signif(c(mu.hat, exp(theta.int)), 5) names(mu) <- c("mu.hat", "mu.L", "mu.U") k.find <- function(p = 0.9, gam = 0.95, sig11 = 1, sig22 = 1, sig12 = 0) { sig1 <- sqrt(sig11) sig2 <- sqrt(sig22) sig1.2 <- sqrt(sig11 - (sig12^2)/sig22) wp <- log(-log(1 - p)) alpha1 <- (1 - gam)/2 beta1 <- gam/2 if (wp >= 0) { v1b <- qnorm(beta1) * sig1 v2b <- qnorm(beta1) * sig2 v1a <- qnorm(1 - alpha1) * sig1 v2a <- qnorm(1 - alpha1) * sig2 } else { v1b <- qnorm(beta1) * sig1 v2b <- qnorm(1 - beta1) * sig2 v1a <- qnorm(1 - alpha1) * sig1 v2a <- qnorm(alpha1) * sig2 } xleft <- v1b - wp * exp(-v2b) + wp xright <- v1a - wp * exp(-v2a) + wp fun0 <- function(xk, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2, gam = gam, wp = wp) { fun <- function(x, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2, gam = gam, wp = wp, xk = xk) { arg <- (xk - wp + wp * exp(-sig2 * x) - (x * sig12)/sig2)/sig1.2 (pnorm(arg) - gam) * dnorm(x) } if (exists("is.R")) { out <- integrate(fun, lower = -6, upper = 6, xk = xk, wp = wp, gam = gam, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2)$value } else { out <- integrate(fun, lower = -6, upper = 6, xk = xk, wp = wp, gam = gam, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2)$integral } out } if (exists("is.R")) { uniroot(fun0, c(xleft, xright), maxiter = 2000, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2, gam = gam, wp = wp)$root } else { uniroot(fun0, c(xleft, xright), sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2, gam = gam, wp = wp)$root } } w.find <- function(hy = 0, gam = 0.95, sig11 = 1, sig22 = 1, sig12 = 0) { sig1 <- sqrt(sig11) sig2 <- sqrt(sig22) sig1.2 <- sqrt(sig11 - (sig12^2)/sig22) alpha1 <- (1 - gam)/2 beta1 <- gam/2 v1a <- qnorm(1 - alpha1) * sig1 v1b <- qnorm(beta1) * sig1 if (v1a + hy >= 0) { v2a <- qnorm(1 - alpha1) * sig2 } else { v2a <- qnorm(alpha1) * sig2 } if (v1b + hy >= 0) { v2b <- qnorm(beta1) * sig2 } else { v2b <- qnorm(1 - beta1) * sig2 } xright <- (v1a + hy) * exp(v2a) xleft <- (v1b + hy) * exp(v2b) fun0 <- function(w, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2, gam = gam, hy = hy) { fun <- function(x, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2, gam = gam, hy = hy, w = w) { arg <- (w * exp(-sig2 * x) - hy - (x * sig12)/sig2)/sig1.2 (pnorm(arg) - gam) * dnorm(x) } if (exists("is.R")) { out <- integrate(fun, lower = -6, upper = 6, w = w, hy = hy, gam = gam, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2)$value } else { out <- integrate(fun, lower = -6, upper = 6, w = w, hy = hy, gam = gam, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2)$integral } out } if (exists("is.R")) { uniroot(fun0, c(xleft, xright), maxiter = 2000, sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2, gam = gam, hy = hy)$root } else { uniroot(fun0, c(xleft, xright), sig2 = sig2, sig12 = sig12, sig1.2 = sig1.2, gam = gam, hy = hy)$root } } if (length(p) > 0) { zp <- rep(z, length(p)) zvecp <- cbind(zp, -zp) xp.hat <- rep(as.numeric(coef.hat %*% cov0), length(p)) + b.hat * wp interval.xp <- xp.hat + zvecp * as.numeric(sqrt(t(cov0) %*% tau.11 %*% cov0 + 2 * t(cov0) %*% tau.12 * wp * b.hat + tau.22 * wp^2 * b.hat^2)) interval.xp <- exp(interval.xp) int.l.mon <- NULL int.u.mon <- NULL for (pj in p) { k.l <- k.find(p = pj, gam = (1 + gam)/2, sig11 = t(cov0) %*% tau.11 %*% cov0/b.hat^2, sig22 = tau.22, sig12 = t(cov0) %*% tau.12/b.hat) k.u <- k.find(p = pj, gam = (1 - gam)/2, sig11 = t(cov0) %*% tau.11 %*% cov0/b.hat^2, sig22 = tau.22, sig12 = t(cov0) %*% tau.12/b.hat) int.l.mon <- c(int.l.mon, exp(u.hat + b.hat * log(-log(1 - pj)) - k.l * b.hat)) int.u.mon <- c(int.u.mon, exp(u.hat + b.hat * log(-log(1 - pj)) - k.u * b.hat)) } xp <- cbind(p, exp(xp.hat), interval.xp, int.l.mon, int.u.mon) dimnames(xp) <- list(NULL, c("p", "xp.hat", "xp.L", "xp.U", "xp.L.m", "xp.U.m")) } if (length(x) > 0) { zp <- rep(z, length(x)) zvecp <- cbind(zp, -zp) x0 <- log(x) phat <- (x0 - u.hat)/b.hat interval.p <- cbind(phat + zp * (1/b.hat) * sqrt(t(cov0) %*% tau.11 %*% cov0 + 2 * t(cov0) %*% tau.12 * (x0 - u.hat) + (x0 - u.hat)^2 * tau.22), phat - zp * (1/b.hat) * sqrt(t(cov0) %*% tau.11 %*% cov0 + 2 * t(cov0) %*% tau.12 * (x0 - u.hat) + (x0 - u.hat)^2 * tau.22)) phat <- 1 - exp(-exp(phat)) interval.p <- 1 - exp(-exp(interval.p)) w.u <- NULL w.l <- NULL sig.u2 <- t(cov0) %*% tau.11 %*% cov0 sig.b2 <- tau.22 * b.hat^2 sig.ub <- t(cov0) %*% tau.12 * b.hat pL.logit <- NULL pU.logit <- NULL for (qi in x) { w.u <- c(w.u, w.find(hy = (log(qi) - u.hat)/b.hat, gam = (1 + gam)/2, sig11 = t(cov0) %*% tau.11 %*% cov0/b.hat^2, sig22 = tau.22, sig12 = t(cov0) %*% tau.12/b.hat)) w.l <- c(w.l, w.find(hy = (log(qi) - u.hat)/b.hat, gam <- (1 - gam)/2, sig11 = t(cov0) %*% tau.11 %*% cov0/b.hat^2, sig22 = tau.22, sig12 = t(cov0) %*% tau.12/b.hat)) zeta <- (log(qi) - u.hat)/b.hat F.hat <- 1 - exp(-exp(zeta)) se.F <- (exp(zeta - exp(zeta))/b.hat) * sqrt(sig.u2 + 2 * zeta * sig.ub + zeta^2 * sig.b2) w <- exp((-z * se.F)/(F.hat * (1 - F.hat))) pL.logit <- c(pL.logit, F.hat/(F.hat + (1 - F.hat) * w)) pU.logit <- c(pU.logit, F.hat/(F.hat + (1 - F.hat)/w)) } tail.prob <- signif(cbind(phat, interval.p, pL.logit, pU.logit, 1 - exp(-exp(w.l)), 1 - exp(-exp(w.u))), 3) tail.prob <- cbind(x, tail.prob) dimnames(tail.prob) <- list(NULL, c("x", "p.hat", "p.L", "p.U", "pL.logit", "pU.logit", "p.L.m", "p.U.m")) } stats <- c(gam, n.time, r) tim <- proc.time()[1:2] - now names(stats) <- c("confidence level", "n observations", "r failures") result <- list(stats = stats, coef = coef, beta = beta, covariate.vector = cov0, mu = mu, xp = xp, tail.prob = tail.prob, time = tim) result } weibull.paper.mon <- function (dat = shockabsorber, gam = 0.95, head = "Weibull Plot", xunit = "cycles", PS = F, tag = "123", stamp.text = " ", choice = 0) { # This function performs a Weibull analysis for right # censored data. # It is assumed that the data frame given to "dat" has two # variables named time and status. status = 1 means that the # corresponding time was the failure time and status = 0 # means that the corresponding time is right censored, i.e., # the actual failure time exceeds the given time. The time # values should all be positive. # gam gives the confidence level for the given bounds, they # are joint confidences for each pair of upper and lower # bounds corresponding to a given value on the probability # scale. # head = "Weibull Plot" modifies the plot title # xunit = "cycles" gives the units for the x axis, # usually "cycles" or "hours" # PS = T means that the plot will be written to # a Postscript file with name: # paste("weibullpaper",tag,sep="") # tag = "123" see PS for usage. # stamp.text = " " is optional text to be put in the lower right # corner of the plot. # choice = 0 indicates that monotone quantile/probability # bounds are desired # 1 indicates that (classical) quantile bounds are # desired # 2 indicates that (classical) probability bounds # are desired # 3 indicates that all three types of bounds are # desired # Fritz Scholz, July 2005 (revised August 2015) check.data <- function (x) { # check.data takes a data frame as input with variables # "time" and "status" # and checks whether the censored data are sufficient for # estimating the Weibull parameters ier <- 1 time1 <- x$time[x$status == 1] time0 <- x$time[x$status == 0] if (length(unique(time1)) > 1) ier <- 0 if (length(unique(time1)) == 1 & length(time0) > 0) { if (max(time0) > max(time1)) ier <- 0 } ier } stamp <- function (msg = date()) { mtext(msg, 1, 4, adj = 1) } ier <- check.data(dat) write(ier, paste("weibull.diagnostics", tag, sep = "")) if (ier == 0) { if (PS == T) { postscript(file = paste("weibullpaper", tag, sep = ""), horizontal = F, height = 8) } prod.lim <- function(times, status) { ord <- order(times) times <- times[ord] status <- status[ord] args <- unique(times[status == 1]) k <- length(args) nn <- length(status) e <- rep(1, nn) cpsurv <- 1 psurv <- 1 cd <- NULL cn <- NULL L.na <- 0 l.na <- 0 for (i in 1:k) { x <- args[i] d <- sum(status[times == x]) n <- sum(e[times >= x]) l.na <- l.na + d/n L.na <- c(L.na, l.na) psurv <- psurv * (1 - d/n) cpsurv <- c(cpsurv, psurv) } L.na <- exp(-L.na) out <- list(args = args, surv.km = cpsurv[-1], surv.na = L.na[-1]) out } time <- dat$time status <- dat$status prob.lines <- c(1e-04, 0.001 * (1:9), 0.01 * (1:9), seq(0.1, 0.9, length.out = 17), 0.9 + 0.01 * (1:9), 0.995, 0.999, 0.9999) l.prob.lines <- log(-log(1 - prob.lines)) prob.label <- c(".001", ".010", ".100", ".200", ".500", ".900", ".990", ".999") y.prob.label <- c(0.001, 0.01, 0.1, 0.2, 0.5, 0.9, 0.99, 0.999) ly.prob.label <- log(-log(1 - y.prob.label)) out1 <- WEIB.base(dat, gam = gam, p = prob.lines, comp.mode = 1) out2 <- prod.lim(time, status) xp.L <- out1$interval.xp[, 2] xp.U <- out1$interval.xp[, 3] out3 <- Weibull.monotone.R(dat, gam = gam, p = prob.lines) xp.L.m <- out3$xp[, 5] xp.U.m <- out3$xp[, 6] mu.hat <- out3$mu[1] mu.L <- out3$mu[2] mu.U <- out3$mu[3] px <- out1$interval.xp[, 1] x <- c(xp.L[xp.L >= 0.5], xp.U, time[status == 1]) x <- sort(x) n <- length(x) epsilon <- 1e-10 l1 <- log10(x[1]) l2 <- log10(x[n]) k1 <- floor(l1 + epsilon) k2 <- ceiling(l2 - epsilon) nlog.scale <- k2 - k1 imax <- 0 if (nlog.scale > 6) { i.count <- 0 n.max.in.range <- 0 while (k1 + i.count + 6 <= k2) { kk1 <- k1 + i.count kk2 <- k1 + i.count + 6 n.in.range <- sum(10^kk1 <= time & time < 10^kk2) if (n.in.range > n.max.in.range) { imax <- i.count n.max.in.range <- n.in.range } i.count <- i.count + 1 } k.temp <- k1 k1 <- k.temp + imax k2 <- k.temp + imax + 6 nlog.scale <- 6 } if (nlog.scale == 1) { m <- floor(k1/3 + epsilon) unit <- floor(10^(3 * m) + epsilon) scale.labels <- as.character(as.integer( floor(10^(k1 - 3 * m) * c(1, 2, 5, 10) + epsilon))) x.scale.labels <- log(floor(10^(k1 - 3 * m) * c(1, 2, 5, 10) + epsilon)) x.lines <- log(10^(k1 - 3 * m) * (1:10)) } if (nlog.scale == 2) { m <- floor((k1 + 1)/3 + epsilon) unit <- floor(10^(3 * m) + epsilon) if (k1 - 3 * m < 0) { scale.labels <- c("0.1", "0.2", "0.5", "1", "2", "5", "10") x.scale.labels <- log(c(0.1, 0.2, 0.5, 1, 2, 5, 10)) } else { scale.labels <- as.character(as.integer(floor( 10^(k1 - 3 * m) * c(1, 2, 5, 10, 20, 50, 100) + epsilon))) x.scale.labels <- log(floor(10^(k1 - 3 * m) * c(1, 2, 5, 10, 20, 50, 100) + epsilon)) } x.lines <- log(10^(k1 - 3 * m) * c(1:9, 10 * (1:10))) } if (nlog.scale == 3) { m <- floor((k1 + 2)/3 + epsilon) unit <- floor(10^(3 * m) + epsilon) if (k1 - 3 * m == -2) { scale.labels <- c("0.01", "0.02", "0.05", "0.1", "0.2", "0.5", "1", "2", "5", "10") x.scale.labels <- log(c(0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10)) } if (k1 - 3 * m == -1) { scale.labels <- c("0.1", "0.2", "0.5", "1", "2", "5", "10", "20", "50", "100") x.scale.labels <- log(c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100)) } if (k1 - 3 * m > -1) { scale.labels <- as.character(as.integer(floor( 10^(k1 - 3 * m) * c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000) + epsilon))) x.scale.labels <- log(floor(10^(k1 - 3 * m) * c(1, 2, 5, 10, 20, 50, 100, 200, 500, 1000) + epsilon)) } x.lines <- log(10^(k1 - 3 * m) * c(1:9, 10 * (1:9), 100 * (1:10))) } if (nlog.scale == 4) { m <- floor((k1 + 3)/3 + epsilon) unit <- floor(10^(3 * m) + epsilon) if (k1 - 3 * m == -3) { scale.labels <- c("0.001", "0.01", "0.1", "1", "10") x.scale.labels <- log(c(0.001, 0.01, 0.1, 1, 10)) } if (k1 - 3 * m == -2) { scale.labels <- c("0.01", "0.1", "1", "10", "100") x.scale.labels <- log(c(0.01, 0.1, 1, 10, 100)) } if (k1 - 3 * m == -1) { scale.labels <- c("0.1", "1", "10", "100", "1000") x.scale.labels <- log(c(0.1, 1, 10, 100, 1000)) } if (k1 - 3 * m > -1) { scale.labels <- as.character(as.integer(floor( 10^(k1 - 3 * m) * c(1, 10, 100, 1000) + epsilon))) x.scale.labels <- log(floor(10^(k1 - 3 * m) * c(1, 10, 100, 1000) + epsilon)) } x.lines <- log(10^(k1 - 3 * m) * c(1:9, 10 * (1:9), 100 * (1:9), 1000 * (1:10))) } if (nlog.scale == 5) { m <- floor((k1 + 3)/3 + epsilon) unit <- floor(10^(3 * m) + epsilon) if (k1 - 3 * m == -3) { scale.labels <- c("0.001", "0.01", "0.1", "1", "10", "100") x.scale.labels <- log(c(0.001, 0.01, 0.1, 1, 10, 100)) } if (k1 - 3 * m == -2) { scale.labels <- c("0.01", "0.1", "1", "10", "100", "1000") x.scale.labels <- log(c(0.01, 0.1, 1, 10, 100, 1000)) } if (k1 - 3 * m == -1) { scale.labels <- c("0.1", "1", "10", "100", "1000", "10000") x.scale.labels <- log(c(0.1, 1, 10, 100, 1000, 10000)) } if (k1 - 3 * m == 0) { scale.labels <- as.character(as.integer(floor( 10^(k1 - 3 * m) * c(1, 10, 100, 1000, 10000) + epsilon))) x.scale.labels <- log(floor(10^(k1 - 3 * m) * c(1, 10, 100, 1000, 10000) + epsilon)) } if (k1 - 3 * m > 0) { scale.labels <- as.character(as.integer(floor( 10^(k1 - 3 * m) * c(1, 10, 100, 1000, 10000, 1e+05) + epsilon))) x.scale.labels <- log(floor(10^(k1 - 3 * m) * c(1, 10, 100, 1000, 10000, 1e+05) + epsilon)) } x.lines <- log(10^(k1 - 3 * m) * c(1:9, 10 * (1:9), 100 * (1:9), 1000 * (1:9), 10000 * (1:10))) } if (nlog.scale == 6) { m <- floor((k1 + 3)/3 + epsilon) unit <- floor(10^(3 * m) + epsilon) if (k1 - 3 * m == -3) { scale.labels <- c("0.001", "0.01", "0.1", "1", "10", "100", "1000") x.scale.labels <- log(c(0.001, 0.01, 0.1, 1, 10, 100, 1000)) } if (k1 - 3 * m == -2) { scale.labels <- c("0.01", "0.1", "1", "10", "100", "1000", "10000") x.scale.labels <- log(c(0.01, 0.1, 1, 10, 100, 1000, 10000)) } if (k1 - 3 * m == -1) { scale.labels <- c("0.1", "1", "10", "100", "1000", "10000", "100000") x.scale.labels <- log(c(0.1, 1, 10, 100, 1000, 10000, 1e+05)) } if (k1 - 3 * m == 0) { scale.labels <- as.character(as.integer(floor( 10^(k1 - 3 * m) * c(1, 10, 100, 1000, 10000, 1e+05) + epsilon))) x.scale.labels <- log(floor(10^(k1 - 3 * m) * c(1, 10, 100, 1000, 10000, 1e+05) + epsilon)) } if (k1 - 3 * m > 0) { scale.labels <- as.character(as.integer(floor( 10^(k1 - 3 * m) * c(1, 10, 100, 1000, 10000, 1e+05, 1e+06) + epsilon))) x.scale.labels <- log(floor(10^(k1 - 3 * m) * c(1, 10, 100, 1000, 10000, 1e+05, 1e+06) + epsilon)) } x.lines <- log(10^(k1 - 3 * m) * c(1:9, 10 * (1:9), 100 * (1:9), 1000 * (1:9), 10000 * (1:9), 1e+05 * (1:10))) } par(pty = "m") if (unit != 1) x.unit <- paste(xunit, " x", as.integer(unit)) else x.unit <- xunit out.of.range <- NULL out.of.range <- c(out.of.range, x[x < 10^k1]) out.of.range <- c(out.of.range, x[x >= 10^k2]) z <- out2$args zz <- z[z >= 10^k1 & z < 10^k2] pz <- 1 - out2$surv.km kpz <- length(pz) if (kpz >= 2) p1 <- c(0, pz[1:(kpz - 1)]) else p1 <- 0 pz <- (pz + p1)/2 lpz <- log(-log(1 - pz)) lzz <- log(zz) lx <- log(x) lppz <- lpz[z >= 10^k1 & z < 10^k2] if (choice == 2 | choice == 3) { low.x <- exp(min(c(min(log(10^k1), lx[1]), max(log(10^k2), lx[n])))) high.x <- exp(max(c(min(log(10^k1), lx[1]), max(log(10^k2), lx[n])))) q0 <- seq(low.x, high.x, (high.x - low.x)/500) out.p <- WEIB.base(dat, gam = gam, p = NULL, q0 = q0, comp.mode = 2, input = out1$out) p.L <- out.p$interval.p[, 2] p.U <- out.p$interval.p[, 3] } plot(c(min(log(10^k1), lx[1]), max(log(10^k2), lx[n])) - log(10^(3 * m)), c(ly.prob.label[1], ly.prob.label[8]), xlab = x.unit, ylab = "probability", axes = F, type = "n") stamp(stamp.text) if (kpz > 30) { points(lzz - log(10^(3 * m)), lppz, col = "purple", pch = 19, cex = 0.3) } else { points(lzz - log(10^(3 * m)), lppz, col = "purple", pch = 19) } axis(1, x.scale.labels, scale.labels, cex = 0.6) axis(2, ly.prob.label, prob.label, srt = 180) mtext(".632", 2, at = 0, cex = 0.6) box() for (i in 1:length(l.prob.lines)) abline(h = l.prob.lines[i]) for (i in 1:length(x.lines)) abline(v = x.lines[i]) abline(h = 0, lty = 2) alpha.mle <- out1$alpha beta.mle <- out1$beta alphaL <- out1$interval.alpha[1] alphaU <- out1$interval.alpha[2] betaL <- out1$interval.beta[1] betaU <- out1$interval.beta[2] lx1 <- min(log(10^k1), lx[1]) - log(10^(3 * m)) + 0.3 * (max(log(10^k2), lx[n]) - min(log(10^k1), lx[1])) ly1 <- ly.prob.label[1] + 0.3 * (ly.prob.label[8] - ly.prob.label[1]) lx2 <- min(log(10^k1), lx[1]) - log(10^(3 * m)) + 0.9 * (max(log(10^k2), lx[n]) - min(log(10^k1), lx[1])) ly2 <- ly.prob.label[1] + 0.2 * (ly.prob.label[8] - ly.prob.label[1]) if (choice == 0) { legend(lx2, ly2, c("m.l.e.", paste(format(signif(100 * gam, 3)), "% monotone qp-confidence bounds")), lty = c(1, 1), xjust = 0.65, yjust = 1, cex = 0.7, col = c("red", "orange"), bg = "white", lwd = 2) } if (choice == 1) { legend(lx2, ly2, c("m.l.e.", paste(format(signif(100 * gam, 3)), "% q-confidence bounds")), lty = c(1, 2), xjust = 0.65, yjust = 1, cex = 0.7, col = c("red", "blue"), bg = "white", lwd = 2) } if (choice == 2) { legend(lx2, ly2, c("m.l.e.", paste(format(signif(100 * gam, 3)), "% p-confidence bounds")), lty = c(1, 4), xjust = 0.65, yjust = 1, cex = 0.7, col = c("red", "darkgreen"), bg = "white", lwd = 2) } if (choice == 3) { legend(lx2, ly2, c("m.l.e.", paste(format(signif(100 * gam, 3)), "% q-confidence bounds"), paste(format(signif(100 * gam, 3)), "% p-confidence bounds"), paste(format(signif(100 * gam, 3)), "% monotone qp-confidence bounds")), lty = c(1, 2, 4, 1), xjust = 0.65, yjust = 1, cex = 0.7, col = c("red", "blue", "darkgreen", "orange"), bg = "white", lwd = 2) } title(head) wx <- (max(log(10^k2), lx[n]) - min(log(10^k1), lx[1])) wy <- ly.prob.label[8] - ly.prob.label[6] rx1 <- min(log(10^k1), lx[1]) - log(10^(3 * m)) rx2 <- rx1 + 0.55 * wx rect(rx1, ly.prob.label[6] - 0.3 * wy, rx2, ly.prob.label[8] + 0.1 * wy, col = "white") text(rx1 + 0.01 * wx, ly.prob.label[8] - 0.12 * wy, expression(hat(alpha)), adj = 0, cex = 0.7, col = "red") text(rx1 + 0.03 * wx, ly.prob.label[8] - 0.18 * wy, paste("= ", format(signif(alpha.mle, 4)), ","), adj = 0, cex = 0.7, col = "red") text(rx1 + 0.14 * wx, ly.prob.label[8] - 0.18 * wy, paste(format(signif(100 * gam, 3)), "% conf. interval (", format(signif(alphaL, 4)), ",", format(signif(alphaU, 4)), ")"), adj = 0, cex = 0.7, col = "blue") text(rx1 + 0.01 * wx, ly.prob.label[8] - 0.35 * wy, expression(hat(beta)), adj = 0, cex = 0.7, col = "red") text(rx1 + 0.03 * wx, ly.prob.label[8] - 0.39 * wy, paste("= ", format(signif(beta.mle, 4)), ","), adj = 0, cex = 0.7, col = "red") text(rx1 + 0.14 * wx, ly.prob.label[8] - 0.39 * wy, paste(format(signif(100 * gam, 3)), "% conf. interval (", format(signif(betaL, 4)), ",", format(signif(betaU, 4)), ")"), adj = 0, cex = 0.7, col = "blue") text(rx1 + 0.01 * wx, ly.prob.label[8] - 0.65 * wy, "MTTF", adj = 0, cex = 0.7, col = "red") text(rx1 + 0.01 * wx, ly.prob.label[8] - 0.85 * wy, expression(hat(mu)), adj = 0, cex = 0.7, col = "red") text(rx1 + 0.03 * wx, ly.prob.label[8] - 0.89 * wy, paste("= ", format(signif(mu.hat, 4)), ","), adj = 0, cex = 0.7, col = "red") text(rx1 + 0.14 * wx, ly.prob.label[8] - 0.89 * wy, paste(format(signif(100 * gam, 3)), "% conf. interval (", format(signif(mu.L, 4)), ",", format(signif(mu.U, 4)), ")"), adj = 0, cex = 0.7, col = "blue") text(rx1 + 0.01 * wx, ly.prob.label[8] - 1.1 * wy, paste("n = ", length(status), ", r = ", sum(status), " failed cases"), adj = 0, cex = 0.7) if (choice == 0) { lines(log(xp.L.m) - log(10^(3 * m)), log(-log(1 - px)), lty = 1, col = "orange", lwd = 2) lines(log(xp.U.m) - log(10^(3 * m)), log(-log(1 - px)), lty = 1, col = "orange", lwd = 2) } if (choice == 1) { lines(log(xp.L) - log(10^(3 * m)), log(-log(1 - px)), lty = 2, col = "blue", lwd = 2) lines(log(xp.U) - log(10^(3 * m)), log(-log(1 - px)), lty = 2, col = "blue", lwd = 2) } if (choice == 2) { lines(log(q0) - log(10^(3 * m)), log(-log(1 - p.L)), lty = 4, col = "darkgreen", lwd = 2) lines(log(q0) - log(10^(3 * m)), log(-log(1 - p.U)), lty = 4, col = "darkgreen", lwd = 2) } if (choice == 3) { lines(log(xp.L) - log(10^(3 * m)), log(-log(1 - px)), lty = 2, col = "blue", lwd = 2) lines(log(xp.U) - log(10^(3 * m)), log(-log(1 - px)), lty = 2, col = "blue", lwd = 2) lines(log(q0) - log(10^(3 * m)), log(-log(1 - p.L)), lty = 4, col = "darkgreen", lwd = 2) lines(log(q0) - log(10^(3 * m)), log(-log(1 - p.U)), lty = 4, col = "darkgreen", lwd = 2) lines(log(xp.L.m) - log(10^(3 * m)), log(-log(1 - px)), lty = 1, col = "orange", lwd = 2) lines(log(xp.U.m) - log(10^(3 * m)), log(-log(1 - px)), lty = 1, col = "orange", lwd = 2) } abline(beta.mle * (log(10^(3 * m)) - log(alpha.mle)), beta.mle, lty = 1, col = "red", lwd = 2) if (PS == T) dev.off() } } WeibullPaper <- function (x = NULL, n = NULL, units = "cycles", LS = T) { # This function gives a Weibull probability plot # for complete or type II censored samples. In the latter # case the vector x contains the lowest r observed failure # times out of n >= r and n needs to be specified. Otherwise # it is assumed that n=r. # The argument LS=T specifies that least squares estimates are # provided and plotted in addition to the maximum likelihood # estimates. # Otherwise no such LS estimates are provided and plotted. # Two types of LS estimates are given, x vs y and y vs x, # as indicated in the list output. # It is assumed that the functions WeibullPaper1, # WeibullPaper2, and WeibullPaper3 are in the work space to # enable 1 cycle log10, 2 cycle log10, and 3 cycle log10 # plotting, chosen to provide a good fit. # If x == NULL the internal default data set below is taken, # see next 3 lines. if (!exists("survreg")) library(survival) if (is.null(x)) x <- c(7, 12.1, 22.8, 23.1, 25.7, 26.7, 29, 29.9, 39.5, 41.9) alpha1 <- NULL beta1 <- NULL alpha2 <- NULL beta2 <- NULL x <- sort(x) r <- length(x) if (is.null(n)) { n <- r } else { if (r > n) { return("must have length of x = r <= n") } } mx <- x[1] Mx <- x[r] Kx <- ceiling(log10(Mx)) kx <- floor(log10(mx)) if (Kx - kx <= 3) unit <- 10^kx if (Kx - kx > 3) unit <- 10^(Kx - 3) p <- (1:r - 0.5)/n wp <- log10(-log(1 - p)) xx <- x/unit if ((Kx - kx) == 1) { if (unit == 1) { units <- units } else { units <- paste(as.integer(unit), units) } WeibullPaper1(units) x0 <- 2.2 L0 <- 10^0.5 } if ((Kx - kx) == 2) { if (unit == 1) { units <- units } else { units <- paste(as.integer(unit), units) } WeibullPaper2(units) x0 <- 1.5 L0 <- 10 } if ((Kx - kx) >= 3) { if (unit == 1) { units <- units } else { units <- paste(as.integer(unit), units) } WeibullPaper3(units) x0 <- 1.5 L0 <- 10^1.5 } x1 <- x0/3 beta.unit <- log10(3) points(xx[xx > 0.8], wp[xx > 0.8], col = "purple", pch = 16, cex = exp(-0.01 * max(n - 20, 0))) if (r < n) { statusx <- c(rep(1, r), rep(0, n - r)) dat.weibull <- data.frame(c(xx, rep(xx[r], n - r)), statusx) } else { statusx <- rep(1, n) dat.weibull <- data.frame(xx, statusx) } names(dat.weibull) <- c("time", "status") out.weibull <- survreg(Surv(time, status) ~ 1, dist = "weibull", data = dat.weibull) alpha.hat <- exp(out.weibull$coef) beta.hat <- 1/out.weibull$scale if (LS == T) { y <- log10(xx) out.py <- lsfit(wp, y) out.yp <- lsfit(y, wp) a.hat <- out.py$coef[1] b.hat <- out.py$coef[2] alpha1 <- 10^a.hat beta1 <- 1/b.hat A.hat <- out.yp$coef[1] B.hat <- out.yp$coef[2] beta2 <- B.hat alpha2 <- 10^(-A.hat/B.hat) } segments(alpha.hat, log10(-log(1 - 1e-04)), alpha.hat, log10(-log(1 - 0.9999)), col = "blue") p1 <- 1e-04 p2 <- 0.9999 wp1 <- log10(-log(1 - p1)) wp2 <- log10(-log(1 - p2)) y1 <- 10^(wp1/beta.hat + log10(alpha.hat)) y2 <- 10^(wp2/beta.hat + log10(alpha.hat)) if (y1 >= 1) { segments(y1, wp1, y2, wp2, col = "red") } else { y1 <- 0 wp1 <- (y1 - log10(alpha.hat)) * beta.hat segments(10^y1, wp1, y2, wp2, col = "red") } if (LS == T) { p1 <- 1e-04 p2 <- 0.9999 wp1 <- log10(-log(1 - p1)) wp2 <- log10(-log(1 - p2)) segments(alpha1, log10(-log(1 - 1e-04)), alpha1, log10(-log(1 - 0.9999)), col = "green") segments(alpha2, log10(-log(1 - 1e-04)), alpha2, log10(-log(1 - 0.9999)), col = "orange") yy1 <- 10^(wp1/beta1 + log10(alpha1)) yy2 <- 10^(wp2/beta1 + log10(alpha1)) if (yy1 >= 1) { segments(yy1, wp1, yy2, wp2, col = "green") } else { yy1 <- 0 wp1 <- (yy1 - log10(alpha1)) * beta1 segments(10^yy1, wp1, yy2, wp2, col = "green") } zz1 <- 10^(wp1/beta2 + log10(alpha2)) zz2 <- 10^(wp2/beta2 + log10(alpha2)) if (zz1 >= 1) { segments(zz1, wp1, zz2, wp2, col = "orange") } else { zz1 <- 0 wp1 <- (zz1 - log10(alpha2)) * beta1 segments(10^zz1, wp1, zz2, wp2, col = "orange") } segments(x0, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995)) - beta1 * beta.unit, col = "green", lty = 2) segments(x0, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995)) - beta2 * beta.unit, col = "orange", lty = 2) } range <- c(mx, Mx) names(range) = c("minimum", "maximum") points(x0, log10(-log(1 - 0.999995)), pch = 16, col = "red") MLE <- c(alpha.hat * unit, beta.hat) names(MLE) <- c("alpha.hat", "beta.hat") LSE.wp.on.log10x <- c(alpha2, beta2) LSE.log10x.on.wp <- c(alpha1, beta1) if (LS == T) { names(LSE.wp.on.log10x) <- c("alpha.hat.2", "beta.hat.2") names(LSE.log10x.on.wp) <- c("alpha.hat.1", "beta.hat.1") } points(x1, log10(-log(1 - 0.999995)) - beta.hat * beta.unit, pch = 16, col = "red") segments(x0, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995)) - beta.hat * beta.unit, col = "red", lty = 2) text(0.7 * L0, log10(-log(1 - 8e-04)), paste("n =", n, ", r =", r), adj = 0) text(0.7 * L0, log10(-log(1 - 0.00045)), substitute(list(hat(alpha)[MLE] == alpha.x, ~~hat(beta)[MLE] == beta.x), list(alpha.x = signif(alpha.hat * unit, 3), beta.x = signif(beta.hat, 3))), adj = 0, col = "red") text(0.7 * L0, log10(-log(1 - 0.00025)), substitute(list(hat(alpha)[LS1] == alpha.x1, ~~hat(beta)[LS1] == beta.x1), list(nx = n, alpha.x1 = signif(alpha1 * unit, 3), beta.x1 = signif(beta1, 3))), adj = 0, col = "green") text(0.7 * L0, log10(-log(1 - 0.00015)), substitute(list(hat(alpha)[LS2] == alpha.x2, ~~hat(beta)[LS2] == beta.x2), list(alpha.x2 = signif(alpha2 * unit, 3), beta.x2 = signif(beta2, 3))), adj = 0, col = "orange") list(range = range, MLE = MLE, LSE.log10x.on.wp = LSE.log10x.on.wp, LSE.wp.on.log10x = LSE.wp.on.log10x) } WeibullPaper1 <- function (units = "cycles") { par(mar = c(1, 1, 1, 0) + 0.1) p <- c(seq(1e-04, 0.001, 1e-04), seq(0.002, 0.01, 0.001), seq(0.02, 0.1, 0.01), seq(0.2, 0.9, 0.1), 0.95, 0.99, 0.995, 0.999, 0.9999) p.marks <- c(1e-04, 5e-04, 0.001, 0.005, 0.01, 0.1, 0.5, 0.9, 0.95, 0.99, 0.999, 0.9999) p.marks.a <- c(".0001", ".0005", ".001", ".005", ".01", ".1", ".5", ".9", ".95", ".99", ".999", ".9999") ticks <- 1:10 tick.marks <- 1:10 tick.marks.a <- as.character(tick.marks) wp.star <- log10(-log(1 - p)) plot(1, 1e-04, axes = F, ylab = "", xlab = units, log = "x", xlim = c(0.7, 10), ylim = c(min(log10(-log(1 - 5e-05)), wp.star), max(wp.star, log10(-log(1 - 0.9999995)))), type = "n") text(0.85, log10(-log(1 - 0.035)), "p", srt = 0, cex = 0.7) text(3.5, log10(-log(1 - 5e-05)), units, cex = 0.9) for (i in 1:length(wp.star)) { segments(1, wp.star[i], 10, wp.star[i], col = "grey") text(0.97, log10(-log(1 - p.marks[i])), p.marks.a[i], srt = 0, cex = 0.6, adj = 1) } for (i in 1:length(ticks)) { segments(ticks[i], log10(-log(0.9999)), ticks[i], log10(-log(1 - 0.9999)), col = "grey") } text(tick.marks, log10(-log(1 - 8e-05)), tick.marks.a, cex = 0.6) x0 <- 2.2 x1 <- x0/3 text(3, log10(-log(1 - 0.99999999)), "Weibull Probability Paper-1 Cycle Log", cex = 1.2) points(x0, log10(-log(1 - 0.999995)), pch = 16) segments(x0, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995))) segments(x1, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995)) - 10 * log10(3)) beta.unit <- log10(3) segments(rep(x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3), rep(0.95 * x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3)) segments(rep(x1, 101), log10(-log(1 - 0.999995)) - (0:100) * log10(3)/10, rep(0.98 * x1, 101), log10(-log(1 - 0.999995)) - (0:100) * log10(3)/10) text(rep(0.945 * x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3), 0:10, cex = 0.5, adj = 1) text(0.9 * x1, log10(-log(1 - 0.999995)) - 0.5 * log10(3), expression(beta), adj = 1, cex = 0.7) segments(1, 0, 10, 0, lty = 2) text(0.97, 0, ".632", adj = 1, cex = 0.6) } WeibullPaper2 <- function (units = "cycles") { par(mar = c(1, 1, 1, 0) + 0.1) p <- c(seq(1e-04, 0.001, 1e-04), seq(0.002, 0.01, 0.001), seq(0.02, 0.1, 0.01), seq(0.2, 0.9, 0.1), 0.95, 0.99, 0.995, 0.999, 0.9999) p.marks <- c(1e-04, 5e-04, 0.001, 0.005, 0.01, 0.1, 0.5, 0.9, 0.95, 0.99, 0.999, 0.9999) p.marks.a <- c(".0001", ".0005", ".001", ".005", ".01", ".1", ".5", ".9", ".95", ".99", ".999", ".9999") ticks <- c(1:9, seq(10, 100, 10)) tick.marks <- c(1:5, seq(10, 50, 10), 100) tick.marks.a <- as.character(tick.marks) wp.star <- log10(-log(1 - p)) plot(1, 1e-04, axes = F, ylab = "", xlab = units, log = "x", xlim = c(0.4, 100), ylim = c(min(log10(-log(1 - 5e-05)), wp.star), max(wp.star, log10(-log(1 - 0.999995)))), type = "n") text(0.7, log10(-log(1 - 0.035)), "p", srt = 0, cex = 0.7) text(10, log10(-log(1 - 5e-05)), units, cex = 0.9) for (i in 1:length(wp.star)) { segments(1, wp.star[i], 100, wp.star[i], col = "grey") text(0.95, log10(-log(1 - p.marks[i])), p.marks.a[i], srt = 0, cex = 0.6, adj = 1) } for (i in 1:length(ticks)) { segments(ticks[i], log10(-log(0.9999)), ticks[i], log10(-log(1 - 0.9999)), col = "grey") } text(tick.marks, log10(-log(1 - 8e-05)), tick.marks.a, cex = 0.6) x0 <- 1.5 x1 <- x0/3 text(10, log10(-log(1 - 0.9999999)), "Weibull Probability Paper-2 Cycle Log", cex = 1.2) points(x0, log10(-log(1 - 0.999995)), pch = 16) segments(x0, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995))) segments(x1, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995)) - 10 * log10(3)) beta.unit <- log10(3) segments(rep(x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3), rep(0.95 * x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3)) segments(rep(x1, 101), log10(-log(1 - 0.999995)) - (0:100) * log10(3)/10, rep(0.98 * x1, 101), log10(-log(1 - 0.999995)) - (0:100) * log10(3)/10) text(rep(0.945 * x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3), 0:10, cex = 0.5, adj = 1) text(0.7 * x1, log10(-log(1 - 0.999995)) - 0.5 * log10(3), expression(beta), adj = 1, cex = 0.7) segments(1, 0, 100, 0, lty = 2) text(0.95, 0, ".632", adj = 1, cex = 0.6) } WeibullPaper3 <- function (units = "cycles") { par(mar = c(1, 1, 1, 0) + 0.1) p <- c(seq(1e-04, 0.001, 1e-04), seq(0.002, 0.01, 0.001), seq(0.02, 0.1, 0.01), seq(0.2, 0.9, 0.1), 0.95, 0.99, 0.995, 0.999, 0.9999) p.marks <- c(1e-04, 5e-04, 0.001, 0.005, 0.01, 0.1, 0.5, 0.9, 0.95, 0.99, 0.999, 0.9999) p.marks.a <- c(".0001", ".0005", ".001", ".005", ".01", ".1", ".5", ".9", ".95", ".99", ".999", ".9999") ticks <- c(1:9, seq(10, 90, 10), seq(100, 1000, 100)) tick.marks <- c(1:5, seq(10, 50, 10), 100, 200, 300, 500, 1000) tick.marks.a <- as.character(tick.marks) wp.star <- log10(-log(1 - p)) plot(1, 1e-04, axes = F, ylab = "", xlab = units, log = "x", xlim = c(0.4, 1000), ylim = c(min(log10(-log(1 - 5e-05)), wp.star), max(wp.star, log10(-log(1 - 0.999995)))), type = "n") text(0.7, log10(-log(1 - 0.035)), "p", srt = 0, cex = 0.7) text(35, log10(-log(1 - 5e-05)), units, cex = 0.9) for (i in 1:length(wp.star)) { segments(1, wp.star[i], 1000, wp.star[i], col = "grey") text(0.95, log10(-log(1 - p.marks[i])), p.marks.a[i], srt = 0, cex = 0.6, adj = 1) } for (i in 1:length(ticks)) { segments(ticks[i], log10(-log(0.9999)), ticks[i], log10(-log(1 - 0.9999)), col = "grey") } text(tick.marks, log10(-log(1 - 8e-05)), tick.marks.a, cex = 0.6) x0 <- 1.5 x1 <- x0/3 text(35, log10(-log(1 - 0.9999999)), "Weibull Probability Paper-3 Cycle Log", cex = 1.2) points(x0, log10(-log(1 - 0.999995)), pch = 16) segments(x0, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995))) segments(x1, log10(-log(1 - 0.999995)), x1, log10(-log(1 - 0.999995)) - 10 * log10(3)) beta.unit <- log10(3) segments(rep(x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3), rep(0.95 * x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3)) segments(rep(x1, 101), log10(-log(1 - 0.999995)) - (0:100) * log10(3)/10, rep(0.98 * x1, 101), log10(-log(1 - 0.999995)) - (0:100) * log10(3)/10) text(rep(0.945 * x1, 11), log10(-log(1 - 0.999995)) - (0:10) * log10(3), 0:10, cex = 0.5, adj = 1) text(0.7 * x1, log10(-log(1 - 0.999995)) - 0.5 * log10(3), expression(beta), adj = 1, cex = 0.7) segments(1, 0, 1000, 0, lty = 2) text(0.95, 0, ".632", adj = 1, cex = 0.6) } WeibullPivot.timing.plot <- function () { # This function produces the plot in Figure 10 of # Inference for the Weibull Distribution. nvec <- c(10, 100, 1000) tn <- c(15.28, 17.78, 56.59) plot(nvec, tn, xlab = "sample size n", ylab = "elapsed time for WeibullPivots(Nsim=10000)", ylim = c(0, 100), cex = 1.6, cex.axis = 1.6, cex.lab = 1.6, pch = 16) ls.out <- lsfit(nvec, tn) abline(ls.out) text(100, 80, "=", adj = 0.5, cex = 1.6) text(100, 70, "=", adj = 0.5, cex = 1.6) text(90, 80, "intercept ", adj = 1, cex = 1.6) text(90, 70, "slope ", adj = 1, cex = 1.6) text(115, 80, format(signif(ls.out$coef[1], 4)), adj = 0, cex = 1.6) text(115, 70, format(signif(ls.out$coef[2], 4)), adj = 0, cex = 1.6) } WeibullPivots <- function(weib.sample = NULL, alpha = 1000, beta = 1.5, n = 10, r = 10, Nsim = 1000, threshold = NULL, graphics = T) { # This function does Weibull analyses for full Weibull samples # of size n, or for type II censored Weibull samples where only # the lowest r sample values are given and the sample size n # and r are specified. # This function WeibullPivots simulates pivot distributions # in the spirit of Bain (1978) "Statistical Analysis of # Reliability and Life-Testing Models", # where tables (based on 10000 simulations) were given for # various pivots, with occasional typos. The accuracy of the # pivot distribution is controlled by Nsim. For Nsim =1000 # and n = 10 it takes about 5 seconds and the computing time # appears to be proportional to Nsim. # Either a Weibull sample is provided as the first argument # weib.sample", or one of size n is generated from # Weibull(scale=alpha, shape=beta). # When taking logarithms of the Weibull data we arrive at a # location scale family of distributions with location parameter # u=log(alpha) and scale b=1/beta. # u.hat and b.hat denote the MLEs of u and b, respectively. # The pivot distributions of u.hat/b.hat and b.hat are shown. # This is followed by a Weibull plot # (based on weibull.paper.mon, # expected to return on exit) with three types of curves, # based on classical methods for quantiles, tail probabilities # and, a newer method that guarantees a monotone set of bounds # in all cases, censored or not. Here it is assumed that the # Weibull sample is uncensored. # Superimposed on this plot are the data points (purple) and # points (black) that represent the 95% confidence intervals for # specific p-quantiles, p= .001, .005, .01, .025, .05, .1, .2 # .3, .4, .5, .6, .7, .8, .95, .975, .99, .995, .999, # based on the Bain pivot method. These bounds have exact # coverage for any n, provided Nsim gets very large. # The other bounds only have approximate coverage, # that gets better as n increases. # Estimates alpha.hat, beta.hat and of the quantiles are # returned, together with corresponding confidence intervals for # various confidence levels, i.e., # gam = .005,.01,.025,.05,.10,.2,.8,.9,.95,.975,.99,.995. # The various confidence levels refer to the individual upper # (lower) bounds", and not to the interval formed by lower and # upper bounds jointly. # Thus 97.5% lower and 97.5% upper confidence bounds L and U # combine to give a 95% confidence interval [L,U] for the # respective target. # Fritz Scholz March 27, 2008 (revised August 5, 2015) if (!exists("survreg")) library(survival) if (is.null(weib.sample) & (alpha <= 0 | beta <= 0)) return("alpha and beta must be >0") if (n - round(n, 0) > 0 | r - round(r, 0) > 0) return("n and r must be integers") if (n < r) return("r cannot be > n") if (r < 2 & is.null(weib.sample)) return("must have r > 1 to have spread in data") if (!is.null(weib.sample)) { if (length(weib.sample) > r) return("length of provided sample cannot exceed r") if (min(weib.sample) <= 0) return("sample data must be positive") if (abs(max(weib.sample) - min(weib.sample)) <= 1e-10) { return("need positive spread in provided sample") } } if (Nsim - round(Nsim, 0) > 0) return(" need integer Nsim") if (Nsim < 1000) return("Nsim >= 1000 recommended (or disable this check)") if (!is.null(threshold)) { if (min(threshold) <= 0) return("need positive thresholds") } gam <- c(0.005, 0.01, 0.025, 0.05, 0.1, 0.2, 0.8, 0.9, 0.95, 0.975, 0.99, 0.995) k <- length(gam)/2 k2 <- length(gam) p = c(0.001, 0.005, 0.01, 0.025, 0.05, seq(0.1, 0.9, 0.1), 0.95, 0.975, 0.99, 0.995, 0.999) statusx <- c(rep(1, r), rep(0, n - r)) if (is.null(weib.sample)) { weibull.sample <- sort(rweibull(n, shape = beta, scale = alpha))[1:r] } else { weibull.sample <- sort(weib.sample)[1:r] } weibull.sample <- c(weibull.sample, rep(weibull.sample[r], n - r)) dat.weibull <- data.frame(weibull.sample, statusx) names(dat.weibull) <- c("time", "status") out.weibull <- survreg(Surv(time, status) ~ 1, dist = "weibull", data = dat.weibull) alpha.hat <- exp(out.weibull$coef) beta.hat <- 1/out.weibull$scale wp <- log(-log(1 - p)) b.hat <- rep(0, Nsim) u.hat <- rep(0, Nsim) for (i in 1:Nsim) { xweib <- sort(rweibull(n, shape = 1, scale = 1)) xweib <- c(xweib[1:r], rep(xweib[r], n - r)) dat <- data.frame(xweib, statusx) out.x <- survreg(Surv(xweib, statusx) ~ 1, dist = "weibull", data = dat) b.hat[i] <- out.x$scale u.hat[i] <- out.x$coef } ub.pivot.quant <- quantile(u.hat/b.hat, gam) b.pivot.quant <- quantile(b.hat, gam) beta.U <- signif(b.pivot.quant * beta.hat, 3) alpha.L <- signif(alpha.hat * exp(-ub.pivot.quant/beta.hat), 5) alpha.hat <- signif(alpha.hat, 5) beta.hat <- signif(beta.hat, 3) xp <- alpha.hat * (-log(1 - p))^(1/beta.hat) xp.L <- NULL for (i in 1:length(p)) { s.gam <- quantile(u.hat/b.hat + log(-log(1 - p[i])) * (1 - 1/b.hat), gam) xp.L <- rbind(xp.L, xp[i] * exp(-s.gam/beta.hat)) } xp.L <- round(xp.L, 1) xp <- round(xp, 1) names(xp) <- paste(p, rep("-quantile", length(p)), sep = "") rownames(xp.L) <- p conf <- rev(rownames(alpha.L)[(k + 1):k2]) alpha1 <- rev(alpha.L[(k + 1):k2]) alpha2 <- alpha.L[1:k] beta1 <- beta.U[1:k] beta2 <- rev(beta.U[(k + 1):k2]) rownames(alpha1) <- conf rownames(alpha2) <- conf alpha.L <- alpha1 alpha.U <- alpha2 beta.L <- beta1 beta.U <- beta2 if (!is.null(threshold)) { TailProb.U <- NULL TailProb.E <- NULL for (i in 1:length(threshold)) { xx <- (log(threshold[i]) - log(alpha.hat)) * beta.hat TailProb.U <- rbind(TailProb.U, quantile(1 - exp(-exp(u.hat + xx * b.hat)), gam)) TailProb.E <- c(TailProb.E, 1 - exp(-(threshold[i]/ alpha.hat)^beta.hat)) } } else { TailProb.E <- NULL TailProb.U <- NULL } if (graphics == T) { par(mfrow = c(2, 1)) hist(u.hat/b.hat, xlab = expression(hat(u)/hat(b)), main = "", nclass = 50, col = c("blue", "orange")) hist(b.hat, xlab = expression(hat(b)), main = "", nclass = 50, col = c("blue", "orange")) readline("hit return\n") par(mfrow = c(1, 1)) unit <- weibull.paper.mon(dat.weibull, choice = 3) points(log(xp.L[, 3]/unit), log(-log(1 - p)), pch = 16) points(log(xp.L[, 10]/unit), log(-log(1 - p)), pch = 16) lines(log(xp.L[, 3]/unit), log(-log(1 - p))) lines(log(xp.L[, 10]/unit), log(-log(1 - p))) points(log(threshold/unit), log(-log(1 - TailProb.U[, 3])), pch = 16, cex = 1.5, col = "red") points(log(threshold/unit), log(-log(1 - TailProb.U[, 10])), pch = 16, cex = 1.5, col = "red") } BOUNDS <- NULL for (i in 1:length(p)) { BOUNDS <- rbind(BOUNDS, rbind(xp.L[i, 12:7], xp.L[i, 1:6])) } TAILP <- NULL for (i in 1:length(threshold)) { TAILP <- rbind(TAILP, rbind(TailProb.U[i, 1:6], TailProb.U[i, 12:7])) } if (!is.null(threshold)) { TAILP <- round(TAILP, 5) TailProb.E <- round(TailProb.E, 5) } BOUNDS <- data.frame(BOUNDS, row.names = paste(as.vector(t(cbind(p, p))), rep(c("-quantile.L", "-quantile.U"), length(p)), sep = "")) if (!is.null(threshold)) { thresh2 <- as.vector(t(cbind(threshold, threshold))) py <- paste(rep("p(", length(thresh2)), thresh2, rep(c(").L ", ").U "), length(threshold)), sep = "") TAILP <- data.frame(TAILP, row.names = py) colnames(TAILP) <- c(paste(100 * gam[12:7], "%", sep = "")) names(TailProb.E) <- paste(rep("p(", length(threshold)), threshold, rep(")", length(threshold)), sep = "") } MLE <- cbind(u.hat, b.hat) colnames(BOUNDS) <- c(paste(100 * gam[12:7], "%", sep = "")) list(alpha.hat = alpha.hat, beta.hat = beta.hat, alpha.beta.bounds = cbind(alpha.L, alpha.U, beta.L, beta.U), p.quantile.estimates = xp, p.quantile.bounds = BOUNDS, Tail.Probability.Estimates = TailProb.E, Tail.Probability.Bounds = TAILP) }