options(warn = -1)
library(fitdistrplus)
## Loading required package: MASS
## Loading required package: survival
library(FAdist)
library(ggplot2)
Sys.setlocale(category = "LC_ALL", locale = "Hebrew")
## [1] "LC_COLLATE=Hebrew_Israel.1255;LC_CTYPE=Hebrew_Israel.1255;LC_MONETARY=Hebrew_Israel.1255;LC_NUMERIC=C;LC_TIME=Hebrew_Israel.1255"
source("build unique WBS code.R")
source("functions.R")
h<-function(time,shape,scale) {
return (dweibull(time,shape,scale)/(1-pweibull(time,shape,scale)) )
}
h3<-function(time,shape,scale,t0) {
return(dweibull3(time,shape,scale,t0)/(1-pweibull3(time,shape,scale,t0)))
}
mean.wei <- function(shape, scale) { return (scale*gamma(1+1/shape) )}
# Source: http://blogs2.datall-analyse.nl/2016/02/17/rcode_three_parameter_weibull/
# Reference: Jerry Lawless, 2003 "Statistical models and methods for lifetime data" (pp. 187-188)
find.t0 <- function(times) {
##determine threshold of 3-parameter Weibull distribution
event<-rep(1,length(times))
#profile function for threshold (=gamma)
profile.t <- function(gammahat, times, event) {
time <- times-gammahat
fit <- survreg(Surv(time, event) ~ 1, dist="weibull", subset=(time>0))
fit$loglik[1]
}
#optimization
#note: for an estimate of gamma which coincides with the maximum in the profile likelihood
#plot (see below) you may sometimes have to vary the value of 1e+6 (e.g., set it to 1e+5)
gammaMLE <- optim(par=0, profile.t, method="L-BFGS-B",
upper=(1-(1/1e+6))*min(times[event==1]),
control=list(fnscale=-1),
times=times, event=event)
gamma <- gammaMLE$par
return(gamma)
}
build.au0<-function() {
a18 <- read.csv("~/1R/Otzma/Otzma 2018 results.csv", stringsAsFactors=FALSE)
a17 <- read.csv("~/1R/Otzma/Otzma 2017 results.csv", stringsAsFactors=FALSE)
a16 <- read.csv("~/1R/Lab/Results.csv", stringsAsFactors = FALSE)
a18$call_date=as.Date.character(a18$call_date,"%d-%m-%y") # only for 2018
a18$start_date=as.Date.character(a18$start_date,"%d-%m-%y") # only for 2018
a18$finish_date=as.Date.character(a18$finish_date,"%d-%m-%y") # only for 2018
#a17$call_date=as.Date.character(a17$call_date) # only for 2017 option 2
#a17$start_date=as.Date.character(a17$start_date) # only for 2017 option 2
#a17$finish_date=as.Date.character(a17$finish_date,"%d/%m/%y") # only for 2017 option 2
a17$call_date=as.Date.character(a17$call_date,"%d-%m-%y") # only for 2017
a17$start_date=as.Date.character(a17$start_date,"%d-%m-%y") # only for 2017
a17$finish_date=as.Date.character(a17$finish_date,"%d-%m-%y") # only for 2017
a16$start_date=as.Date.numeric(a16$Event_Date, origin = "1970-01-01") # only for history results
au0 <- a17[,c(3,5:19,21)] # for 2017
au0 = rbind(au0,a18[,c(3,5:19,21)]) # add 2018
return(au0)
}
hz.plot<-function(times) {
n=length(times)
#print(n)
Prob=c(0.5*((1:n)/n+(0:(n-1))/n))
yi=log(log((1-Prob)^-1))
par(yaxt="n")
plot(x=sort(times), y=yi, log="x",xlim=range(times), type='b', lwd=2,
main=paste("Failure Rate Distribution"),
ylab="Failure Cumulative Distribution [%]",
xlab="Rate [Days]")
par(yaxt="s")
seqn = 1
if (length(yi)>11) seqn = 2
if (length(yi)>17) seqn = 3
axis(side = 2, at = yi[seq(1,length(yi),seqn)],labels = format(Prob[seq(1,length(yi),seqn)]*100,digits = 2),tcl=0,tck=1, lty = 'dashed',lwd = 0.2, cex.axis=(0.5+1/length(yi)),las=1)
legend(x="topleft",legend = c("Empirical","Weibull 2-par fit","Weibull 3-par fit")
, col=c('black','blue','green'),lty=c(1,2,4))
t0=find.t0(times)
names(t0)="t0"
q3=fitdist(data = c(times-t0+1), distr = "weibull")
q2=fitdist(data = times ,distr = "weibull")
fit=q2$estimate
print(fit)
sd=q2$sd
curve(expr = log(log((1-pweibull(x,fit[1],fit[2]))^-1)), xlim = c(range(times)), add = T, col="blue", lwd=2, lty = 2)
curve(expr = log(log((1-pweibull(x,fit[1]+sd[1],fit[2]+sd[2]))^-1)), xlim = c(range(times)), add = T, col="blue", lwd=0.5, lty=2)
curve(expr = log(log((1-pweibull(x,fit[1]-sd[1],fit[2]-sd[2]))^-1)), xlim = c(range(times)), add = T, col="blue", lwd=0.5, lty=2)
hz2=h(times,fit[1],fit[2])
k2=ks.test(Prob,pweibull(times,fit[1],fit[2]),exact = F)$p.value
names(k2)="k2.p.val"
print(k2)
fit=q3$estimate
sd=q3$sd
curve(expr = log(log((1-pweibull3(x,fit[1],fit[2],t0))^-1)), xlim = range(times-t0+1), add = T, col="green", lwd=2, lt = 4)
curve(expr = log(log((1-pweibull3(x,fit[1]+sd[1],fit[2]+sd[2],t0))^-1)), xlim = range(times-t0+1), add = T, col="green", lwd=0.5, lty=4)
curve(expr = log(log((1-pweibull3(x,fit[1]-sd[1],fit[2]-sd[2],t0))^-1)), xlim = range(times-t0+1), add = T, col="green", lwd=0.5, lty=4)
hz3=h3(times,fit[1],fit[2],t0)
print(fit)
print(ceiling(t0))
k3=ks.test(Prob,pweibull3(times,fit[1],fit[2],t0),exact = F)$p.value
names(k3)="k3.p.val"
print(k3)
d=data.frame(time=times,hz2=hz2,hz3=hz3,shape2=q2$estimate[1],scale2=q2$estimate[2],shape.sd=q2$sd[1],scale.sd=q2$sd[2],
shape3=fit[1],scale3=fit[2],shape3.sd=sd[1],scale3.sd=sd[2],t0,ll2=q2$loglik,ll3=q3$loglik,k2,k3) #shape=fit[1],scale=fit[2],shape.sd=q$sd[1],scale.sd=q$sd[2],
return(d)
}
au0=build.au0()
l=build.WBS.code(au0 = au0)
code=l$code
au0=l$au0
df=NULL
#c=code$code[29]
for (c in levels(code$code)) {
ix=which(code$code == c)
print(paste("code ix:",ix))
print(paste("code:",c))
tz=code$origin[ix]
if (is.na(tz) | tz==0) next
au2=au0[au0$code==c,]
au2=au2[order(au2$start_date),]
if (length(au2$machine)<4) next
m=build_dt1( au2 )
#m = m[order(m$start_date),]
n=length(m$dt)
range=c(floor(as.numeric(min(m$start_date))/365)+1970, ceiling(as.numeric(max(m$start_date))/365)+1970)
range=paste0(range,"-01-01")
range=as.Date.character(range)
rate.start=2
if (tz>=as.numeric(range[1])) rate.start=3
par(xaxt="n")
plot(x=m$start_date[rate.start:n], y=m$dt[rate.start:n], xlim=range,
type='b',main=paste("Failure rate vs Time for",c),
ylab="Rate [days]" ,xlab="")
par(xaxt="s")
seqn = 1
if (n>5) seqn = 2
if (n>11) seqn = 3
axis(side = 1,at = m$start_date[seq(1,n,seqn)], tick = T,
labels = as.character(m$start_date[seq(1,n,seqn)],format("%d-%m")),
las=3, cex.axis = (0.5+1/length(m$start_date)), hadj = 1.7)
m$start_date= (as.numeric(m$start_date)) #-tz
fr = data.frame(time = (m$start_date) , rate = m$dt)
fr=fr[rate.start:n,]
d = data.frame(ix=ix,code=as.character(c),hz.plot(fr$rate))
df= rbind(df,d)
plot.new()
}
## [1] "code ix: 2"
## [1] "code: "
## [1] "code ix: 17"
## [1] "code: auto0מכונות"
## [1] "code ix: 14"
## [1] "code: auto1בקרות"
## shape scale
## 0.8504077 57.6905762
## k2.p.val
## 0.9962552
## shape scale
## 0.5849338 38.5913947
## t0
## 7
## k3.p.val
## 0.9962552
## [1] "code ix: 8"
## [1] "code: auto1חשמל"
## shape scale
## 1.436167 31.948633
## k2.p.val
## 0.9996333
## shape scale
## 1.277355 28.987439
## t0
## 3
## k3.p.val
## 0.9996333
## [1] "code ix: 1"
## [1] "code: auto1מכונות"
## shape scale
## 0.9921532 8.6983956
## k2.p.val
## 0.617195
## shape scale
## 0.9921535 8.6983972
## t0
## 1
## k3.p.val
## 0.2301328
## [1] "code ix: 16"
## [1] "code: auto1מסגרות"
## [1] "code ix: 7"
## [1] "code: auto2בקרות"
## shape scale
## 0.5363409 58.5198257
## k2.p.val
## 0.9996333
## shape scale
## 0.5363409 58.5198312
## t0
## 1
## k3.p.val
## 0.9996333
## [1] "code ix: 15"
## [1] "code: auto2חשמל"
## shape scale
## 0.8662477 52.5471828
## k2.p.val
## 0.8927783
## shape scale
## 0.5283845 22.4357288
## t0
## 14
## k3.p.val
## 0.8927783
## [1] "code ix: 5"
## [1] "code: auto2מכונות"
## shape scale
## 1.092037 20.605040
## k2.p.val
## 0.9999968
## shape scale
## 1.092037 20.605042
## t0
## 1
## k3.p.val
## 0.9919639
## [1] "code ix: 10"
## [1] "code: auto2מסגרות"
## shape scale
## 1.978005 63.379084
## k2.p.val
## 1
## shape scale
## 1.149547 44.852493
## t0
## 14
## k3.p.val
## 1
## [1] "code ix: 3"
## [1] "code: auto3בקרות"
## shape scale
## 0.5427151 33.2745163
## k2.p.val
## 0.9639452
## shape scale
## 0.5427152 33.2745201
## t0
## 1
## k3.p.val
## 0.9639452
## [1] "code ix: 6"
## [1] "code: auto3חשמל"
## shape scale
## 0.950746 43.239959
## k2.p.val
## 0.993356
## shape scale
## 0.7835406 36.1048142
## t0
## 4
## k3.p.val
## 0.993356
## [1] "code ix: 45"
## [1] "code: auto3כללי"
## [1] "code ix: 4"
## [1] "code: auto3מכונות"
## shape scale
## 1.04472 17.13479
## k2.p.val
## 0.9951661
## shape scale
## 1.044721 17.134796
## t0
## 1
## k3.p.val
## 0.9180184
## [1] "code ix: 18"
## [1] "code: auto3מסגרות"
## shape scale
## 0.4367336 60.5920570
## k2.p.val
## 0.9639452
## shape scale
## 0.4367337 60.5920671
## t0
## 1
## k3.p.val
## 0.9639452
## [1] "code ix: 9"
## [1] "code: auto4בקרות"
## shape scale
## 0.7747975 46.0443660
## k2.p.val
## 0.9793631
## shape scale
## 0.7747976 46.0443687
## t0
## 1
## k3.p.val
## 0.9793631
## [1] "code ix: 11"
## [1] "code: auto4חשמל"
## shape scale
## 1.08337 45.04122
## k2.p.val
## 0.7590978
## shape scale
## 1.08337 45.04122
## t0
## 1
## k3.p.val
## 0.7590978
## [1] "code ix: 13"
## [1] "code: auto4מכונות"
## shape scale
## 0.8634671 32.6702019
## k2.p.val
## 0.9978966
## shape scale
## 0.8634673 32.6702035
## t0
## 1
## k3.p.val
## 0.9978966
## [1] "code ix: 12"
## [1] "code: auto4מסגרות"
## shape scale
## 11.29599 121.39137
## k2.p.val
## 0.9962552
## shape scale
## 8143.076 79362.830
## t0
## -79239
## k3.p.val
## 0.9962552
## [1] "code ix: 35"
## [1] "code: bullmer0בקרות"
## shape scale
## 1.085926 40.457380
## k2.p.val
## 1
## shape scale
## 0.8390584 32.1964668
## t0
## 5
## k3.p.val
## 0.9639452
## [1] "code ix: 34"
## [1] "code: bullmer0חשמל"
## shape scale
## 0.899048 23.030611
## k2.p.val
## 0.9999652
## shape scale
## 0.8990482 23.0306130
## t0
## 1
## k3.p.val
## 0.9780359
## [1] "code ix: 36"
## [1] "code: bullmer0מכונות"
## shape scale
## 0.7742318 37.6910408
## k2.p.val
## 0.9882611
## shape scale
## 0.774232 37.691043
## t0
## 1
## k3.p.val
## 0.9882611
## [1] "code ix: 44"
## [1] "code: c-scan1בקרות"
## [1] "code ix: 47"
## [1] "code: c-scan1חשמל"
## shape scale
## 2.195374 21.867571
## k2.p.val
## 0.9962552
## shape scale
## 1.032311 13.476935
## t0
## 7
## k3.p.val
## 0.9962552
## [1] "code ix: 48"
## [1] "code: c-scan1כללי"
## [1] "code ix: 43"
## [1] "code: c-scan1מכונות"
## shape scale
## 1.34744 45.25426
## k2.p.val
## 1
## shape scale
## 0.9322998 35.6593944
## t0
## 6
## k3.p.val
## 0.9882611
## [1] "code ix: 42"
## [1] "code: c-scan1מסגרות"
## shape scale
## 0.6247009 44.7049168
## k2.p.val
## 0.9375027
## shape scale
## 0.624701 44.704919
## t0
## 1
## k3.p.val
## 0.9375027
## [1] "code ix: 21"
## [1] "code: eagle0בקרות"
## shape scale
## 0.8191674 20.4675735
## k2.p.val
## 0.9639452
## shape scale
## 0.6849288 15.8870831
## t0
## 3
## k3.p.val
## 0.9639452
## [1] "code ix: 19"
## [1] "code: eagle0חשמל"
## shape scale
## 4.219727 118.237821
## k2.p.val
## 0.9962552
## shape scale
## 0.74634 25.26493
## t0
## 79
## k3.p.val
## 0.9962552
## [1] "code ix: 20"
## [1] "code: eagle0מכונות"
## shape scale
## 0.877023 14.347358
## k2.p.val
## 0.958596
## shape scale
## 0.8770231 14.3473599
## t0
## 1
## k3.p.val
## 0.6070057
## [1] "code ix: 31"
## [1] "code: eastman0בקרות"
## shape scale
## 1.003677 14.023454
## k2.p.val
## 0.674894
## shape scale
## 0.8851663 12.1756928
## t0
## 2
## k3.p.val
## 0.674894
## [1] "code ix: 33"
## [1] "code: eastman0חשמל"
## [1] "code ix: 32"
## [1] "code: eastman0מכונות"
## shape scale
## 0.7026433 24.1444340
## k2.p.val
## 0.9793631
## shape scale
## 0.7026434 24.1444368
## t0
## 1
## k3.p.val
## 0.9793631
## [1] "code ix: 41"
## [1] "code: matec0בקרות"
## shape scale
## 1.31296 26.51334
## k2.p.val
## 1
## shape scale
## 1.090065 23.128864
## t0
## 3
## k3.p.val
## 1
## [1] "code ix: 38"
## [1] "code: matec0חשמל"
## shape scale
## 1.04520 44.81694
## k2.p.val
## 1
## shape scale
## 1.045201 44.816943
## t0
## 1
## k3.p.val
## 1
## [1] "code ix: 37"
## [1] "code: matec0מכונות"
## shape scale
## 1.043526 39.983139
## k2.p.val
## 0.9639452
## shape scale
## 1.043526 39.983140
## t0
## 1
## k3.p.val
## 0.9639452
## [1] "code ix: 40"
## [1] "code: matec0מסגרות"
## [1] "code ix: 39"
## [1] "code: matec0צנרת"
## [1] "code ix: 26"
## [1] "code: matrix1בקרות"
## shape scale
## 1.100579 138.911237
## k2.p.val
## 0.9962552
## shape scale
## 0.6092743 99.2225501
## t0
## 11
## k3.p.val
## 0.9962552
## [1] "code ix: 24"
## [1] "code: matrix1חשמל"
## shape scale
## 0.9679526 56.7467799
## k2.p.val
## 0.9375027
## shape scale
## 0.7633709 46.4568621
## t0
## 5
## k3.p.val
## 0.9999997
## [1] "code ix: 25"
## [1] "code: matrix1כללי"
## shape scale
## 1.334752 70.536945
## k2.p.val
## 0.9999968
## shape scale
## 1.003147 61.579260
## t0
## 5
## k3.p.val
## 0.8927783
## [1] "code ix: 22"
## [1] "code: matrix1מכונות"
## shape scale
## 1.032531 9.528735
## k2.p.val
## 0.6993742
## shape scale
## 1.032532 9.528736
## t0
## 1
## k3.p.val
## 0.3802429
## [1] "code ix: 23"
## [1] "code: matrix1מסגרות"
## shape scale
## 0.8363598 127.3772728
## k2.p.val
## 0.9962552
## shape scale
## 0.5948187 102.7045033
## t0
## 5
## k3.p.val
## 0.9962552
## [1] "code ix: 46"
## [1] "code: matrix2בקרות"
## [1] "code ix: 30"
## [1] "code: matrix2חשמל"
## shape scale
## 0.6722939 30.2753450
## k2.p.val
## 0.9962552
## shape scale
## 0.6722941 30.2753480
## t0
## 1
## k3.p.val
## 0.9962552
## [1] "code ix: 27"
## [1] "code: matrix2כללי"
## shape scale
## 3.09157 92.87408
## k2.p.val
## 0.9639452
## shape scale
## 0.5813966 21.8541182
## t0
## 52
## k3.p.val
## 0.9639452
## [1] "code ix: 29"
## [1] "code: matrix2מכונות"
## shape scale
## 0.9796922 8.1752733
## k2.p.val
## 0.8775343
## shape scale
## 0.9796925 8.1752749
## t0
## 1
## k3.p.val
## 0.6487256
## [1] "code ix: 28"
## [1] "code: matrix2מסגרות"
## shape scale
## 1.131424 22.532695
## k2.p.val
## 0.9999652
## shape scale
## 0.8524051 18.2638580
## t0
## 3
## k3.p.val
## 0.8186212
f = file("fr vs time v2 2-par.csv","wt")
write.csv(df,f)
close(f)