# It is easiest to think about residual life as starting at the origin,
# but for this example, let us set the lower limit at 100 days.
A <- vec2par(c(100, 2649, 2.11), type="gov")
f <- 0.47 # Both computations of Phi show 0.6455061
"afunc" <- function(p) { return(par2qua(p,A,paracheck=FALSE)) }
tmpa <- 1/cmlmomco(f=0, A); tmpb <- (1-f)*par2qua(f,A,paracheck=FALSE)
Phiu1 <- tmpa * ( tmpb + integrate(afunc,0,f)$value )
Phiu2 <- stttlmomco(f, A)
if (FALSE) {
# The TTT-plot (see Nair et al. (2013, p. 173))
n <- 30; X <- sort(rlmomco(n, A)); lmr <- lmoms(X) # simulated lives and their L-moments
# recognize here that the "fit" is to the lifetime data themselves and not to special
# curves or projections of the data to other scales
"Phir" <- function(r, X, sort=TRUE) {
n <- length(X); if(sort) X <- sort(X)
if(r == 0) return(0) # can use 2:r as X_{0:n} is zero
Tau.rOFn <- sapply(1:r, function(j) { Xlo <- ifelse((j-1) == 0, 0, X[(j-1)]);
return((n-j+1)*(X[j] - Xlo)) })
return(sum(Tau.rOFn))
}
Xbar <- mean(X); rOFn <- (1:n)/n # Nair et al. (2013) are clear r/n used in the Phi(u)
Phi <- sapply(1:n, function(r) { return(Phir(r,X, sort=FALSE)) }) / (n*Xbar)
layout(matrix(1:3, ncol=1))
plot(rOFn, Phi, type="b",
xlab="NONEXCEEDANCE PROBABILITY", ylab="SCALED TOTAL TIME ON TEST")
lines(rOFn, stttlmomco(rOFn, A), lwd=2, col=8) # solid grey, the parent distribution
par1 <- pargov(lmr); par2 <- pargov(lmr, xi=min(X)) # notice attempt to "fit at minimum"
lines(pp(X), stttlmomco(rOFn, par1)) # now Weibull (i/(n+1)) being used for F via pp()
lines(pp(X), stttlmomco(rOFn, par2), lty=2) # perhaps better, but could miss short lives
F <- nonexceeds(f01=TRUE)
plot(pp(X), sort(X), xlab="NONEXCEEDANCE PROBABILITY", ylab="TOTAL TIME ON TEST (DAYS)")
lines(F, qlmomco(F, A), lwd=2, col=8) # the parent again
lines(F, qlmomco(F, par1), lty=1); lines(F, qlmomco(F, par2), lty=2) # two estimated fits
plot(F, lrzlmomco(F, par2), col=2, type="l") # Lorenz curve from L-moment fit (red)
lines(F, bfrlmomco(F, par2), col=3, lty=2) # Bonferroni curve from L-moment fit (green)
lines(F, lkhlmomco(F, par2), col=4, lty=4) # Leimkuhler curve from L-moment fit (blue)
lines(rOFn, Phi) # Scaled Total Time on Test
}
Run the code above in your browser using DataLab