R/BBB.r

In WeibullR: Weibull Analysis for Reliability Engineering

BBB<-function(xdata, xfit,n=NULL, CI=.90, unrel=NULL, type="horizontal", nknots=NULL)  {						
##  xdata is  a list containing either dpoints and/or dlines such as available from a wblr object						
## xfit is a list containing dist and fit elements, both can be extracted from a wblr.fit object, or generated by preprocessing.						
						
	if(is(xdata,"list"))  {					
	## validate dpoints and/or dlines elements by comparing column names					
		if(is.null(xdata$dpoints) && is.null(xdata$dlines))  {				
			stop("xdata does not contain dpoints or dlines elements")			
		}else{				
			if(!is.null(xdata$dpoints))  {			
				if(!all.equal(names(xdata$dpoints), c("time","ppp","adj_rank","weight"))) {		
					stop("xdata$dpoints improperly formed")	
				}		
			}			
			if(!is.null(xdata$dlines))  {			
				if(!all.equal(names(xdata$dpoints), c("t1", "t2", "ppp","adj_rank","weight")))  {		
					stop("xdata$dlines improperly formed")	
				}		
			}			
			if(!is.null(xdata$lrq_frame))  {			
				if(is.null(n))  n <- sum(xdata$lrq_frame$qty)		
			}			
		}				
	}else{					
		stop("xdata must be a list such as available from a wblr object")				
	}					
						
	if(is(xfit,"list"))  {					
		if(is.null(xfit$dist) || is.null(xfit$fit))  {				
			stop("xfit must contain both  dist and fit objects")			
		}else{				
		## It turns out that this code is similar to all fitting methods:				
			#dist<-xfit$dist			
			if(any(tolower(xfit$dist) %in% c("weibull","weibull2p","weibull3p","weibull3"))){			
				base_distribution<-"weibull"		
				npar<-2					
				if(any(tolower(xfit$dist) %in% c("weibull3p","weibull3"))){		
					base_distribution<-"weibull"	
					npar<-3	
				}		
			}else{			
				if(any(tolower(xfit$dist) %in% c("lnorm", "lognormal","lognormal2p", "lognormal3p", "lognormal3"))){		
					base_distribution<-"lnorm"	
					npar<-2	
					if(any(tolower(xfit$dist) %in% c( "lognormal3p", "lognormal3"))){	
						base_distribution<-"lognormal"
						npar<-3
					}	
						
				}else{		
		#			if(!xfit$dist=="gumbel") {	
				## Note:  only lslr contains experimental support for "gumbel"		
					stop(paste0("dist argument ", xfit$dist, "is not recognized"))	
		#			}	
				}		
						
			}			
## Perhaps a test that the xfit$fit is indeed a numeric vector would be better						
			fit_par<-xfit$fit			
			if(!is.null(xfit$n))  {			
				if(!is.null(attr(fit_par,"data_types")) ) {		
					n<-attr(fit_par,"data_types")[1]	
				}		
			}			
		}				
	}					
	if(is.null(n))  {					
		stop("must have n as total failures and suspensions")				
	}					
						
	if(length(unrel)>0)  {					
	dp<-unrel					
	}else{					
	## these descriptive percentiles match Minitab unchangeable defaults					
	dp=c(seq(.01,.09,by=.01),seq(.10,.90,by=.10),seq(.91,.99, by=.01))					
	}					
						
## Need to combine ppp and adjusted ranks for points and lines.						
	sx<-NULL					
	if(!is.null(xdata$dpoints)) {					
	sx<-xdata$dpoints[,2:4]					
	}					
	if(!is.null(xdata$dlines)) {					
	sx<-rbind(sx,xdata$dlines[,3:5])					
	sx<-sx[order(sx$adj_rank),]					
	}					
						
## Beta Binomial "Z" factors are non-parametric						
	Zlo<-qbeta((1-CI)/2,sx$adj_rank,n-sx$adj_rank+1)					
	Zhi<-qbeta(1-(1-CI)/2,sx$adj_rank,n-sx$adj_rank+1)					
						
## identify the bound limit points  and add them to the descriptive percentiles						
	h.min<-min(sx$ppp)					
	h.max<-max(sx$ppp)					
	v.min.lo<-min(Zhi)					
	v.max.lo<-max(Zhi)					
	v.min.up<-min(Zlo)					
	v.max.up<-max(Zlo)					
						
	dp<-c(dp, h.min, h.max, v.min.lo, v.max.lo, v.min.up, v.max.up)					
	dp <- unique(signif(dp[order(dp)]))					
# signif() has been used to eliminate any identical looking descriptive						
# percentiles that differ only at place far from the decimal point						
						
## get the position in the dp vector for the bound limit points						
	h.min.pos<-match(signif(h.min), dp)
	h.max.pos<-match(signif(h.max), dp)
	v.min.lo.pos<-match(signif(v.min.lo), dp)
	v.max.lo.pos<-match(signif(v.max.lo), dp)
	v.min.up.pos<-match(signif(v.min.up), dp)
	v.max.up.pos<-match(signif(v.max.up), dp)
					
						
	if(npar==2)  {					
		if(base_distribution=="weibull")  {				
			lower.h<- data.frame(p=sx$ppp,lower=qweibull(Zlo,fit_par[2],fit_par[1]))			
			lower.v<- data.frame(p=Zhi, lower=qweibull(sx$ppp, fit_par[2],fit_par[1]))			
			fit.dp<- qweibull(dp,fit_par[2],fit_par[1])			
			upper.h<- data.frame(p=sx$ppp, upper=qweibull(Zhi,fit_par[2],fit_par[1]))			
			upper.v<-data.frame(p=Zlo, upper=qweibull(sx$ppp, fit_par[2],fit_par[1]))			
		}else{				
			if(base_distribution=="lognormal")  {			
				lower.h<- data.frame(p=sx$ppp,lower=qlnorm(Zlo,fit_par[1],fit_par[2]))		
				lower.v<- data.frame(p=Zhi, lower=qlnorm(sx$ppp, fit_par[1],fit_par[2]))		
				fit.dp<- qlnorm(dp,fit_par[1],fit_par[2])		
				upper.h<- data.frame(p=sx$ppp, upper=qlnorm(Zhi,fit_par[1],fit_par[2]))		
				upper.v<-data.frame(p=Zlo, upper=qlnorm(Zhi,fit_par[1],fit_par[2]))		
			}else{			
				stop(paste0("distribution ", base_distribution, " not supported for bbb."))		
			}			
						
		}				
	}					
						
	if(npar==3)  {					
		if(base_distribution=="weibull")  {				
			lower.h<- data.frame(p=sx$ppp,lower=qweibull(Zlo,fit_par[2],fit_par[1]) + fit_par[3])			
			lower.v<- data.frame(p=Zhi, lower=qweibull(sx$ppp, fit_par[2],fit_par[1]) + fit_par[3])			
			fit.dp<- qweibull(dp,fit_par[2],fit_par[1]) + fit_par[3]			
			upper.h<- data.frame(p=sx$ppp, upper=qweibull(Zhi,fit_par[2],fit_par[1]) + fit_par[3])			
			upper.v<-data.frame(p=Zlo, upper=qweibull(sx$ppp, fit_par[2],fit_par[1]) + fit_par[3])			
		}else{				
			if(base_distribution=="lognormal")  {			
				lower.h<- data.frame(p=sx$ppp,lower=qlnorm(Zlo,fit_par[1],fit_par[2]) + fit_par[3])	
				lower.v<- data.frame(p=Zhi, lower=qlnorm(sx$ppp, fit_par[1],fit_par[2]) + fit_par[3])		
				fit.dp<- qlnorm(dp,fit_par[1],fit_par[2]) + fit_par[3]		
				upper.h<- data.frame(p=sx$ppp, upper=qlnorm(Zhi,fit_par[1],fit_par[2]) + fit_par[3])		
				upper.v<-data.frame(p=Zlo, upper=qlnorm(Zhi,fit_par[1],fit_par[2]) + fit_par[3])		
			}else{			
				stop(paste0("distribution ", base_distribution, " not supported for bbb."))		
			}			
		}				
	}					
						
if(is.null(nknots))  {						
	if(type=="horizontal")  {					
		lo.spline<-smooth.spline(lower.h)				
		up.spline<-smooth.spline(upper.h)				
	}					
						
	if(type=="vertical")  {					
		lo.spline<-smooth.spline(lower.v)				
		up.spline<-smooth.spline(upper.v)				
	}					
						
						
	if(type %in% c("valid", "extrapolated"))  {					
## prepare combined h and v prep points for both lower and upper						
	lower.hv<-rbind(lower.h, lower.v)					
	lower.hv<-lower.hv[order(lower.hv$p),]					
	upper.hv<-rbind(upper.h, upper.v)					
	upper.hv<-upper.hv[order(upper.hv$p),]					
						
	lo.spline<-smooth.spline(lower.hv)					
	up.spline<-smooth.spline(upper.hv)					
	}					
}else{
	if(type=="horizontal")  {					
		lo.spline<-smooth.spline(lower.h, nknots=nknots)				
		up.spline<-smooth.spline(upper.h, nknots=nknots)				
	}					
						
	if(type=="vertical")  {					
		lo.spline<-smooth.spline(lower.v, nknots=nknots)				
		up.spline<-smooth.spline(upper.v, nknots=nknots)				
	}					
						
						
	if(type %in% c("valid", "extrapolated"))  {					
## prepare combined h and v prep points for both lower and upper						
	lower.hv<-rbind(lower.h, lower.v)					
	lower.hv<-lower.hv[order(lower.hv$p),]					
	upper.hv<-rbind(upper.h, upper.v)					
	upper.hv<-upper.hv[order(upper.hv$p),]					
						
	lo.spline<-smooth.spline(lower.hv, nknots=nknots)					
	up.spline<-smooth.spline(upper.hv, nknots=nknots)
	}
}						
	lo.list<- predict(lo.spline,dp)					
	up.list<- predict(up.spline, dp)					
						
	bounds<-data.frame(Prob=dp, Lower=lo.list$y, Fit=fit.dp, Upper=up.list$y)					
						
## now eliminate invalid points from the extrapolated bounds						
						
	if(type == "horizontal")  {					
		if(h.min.pos > 1)  {				
			bounds$Lower[1:(h.min.pos-1)] <- NA			
			bounds$Upper[1:(h.min.pos-1)] <- NA			
		}				
		if(h.max.pos < length(dp) )  {				
			bounds$Lower[(h.max.pos+1): length(dp)] <- NA			
			bounds$Upper[(h.max.pos+1): length(dp)] <- NA			
		}				
	}					
						
	if(type == "vertical")  {					
		bounds$Lower[1:(v.min.lo.pos-1)] <- NA				
		if(v.max.lo.pos < length(dp) )  bounds$Lower[(v.max.lo.pos+1):length(dp)] <- NA				
		if(v.min.up.pos >1 )  bounds$Upper[1:(v.min.up.pos-1)] <- NA				
		bounds$Upper[(v.max.up.pos+1): length(dp)] <- NA				
	}					
						
	if(type == "valid")  {					
		if(h.min.pos > 1)  bounds$Lower[1:(h.min.pos-1)] <- NA				
		if(v.max.lo.pos < length(dp) )  bounds$Lower[(v.max.lo.pos+1): length(dp)] <- NA				
		if(v.min.up.pos >1 )  bounds$Upper[1:(v.min.up.pos-1)] <- NA				
		if(h.max.pos < length(dp) )  bounds$Upper[(h.max.pos+1): length(dp)] <- NA				
	}					
						
						
	return(bounds)					
						
}