library(ggplot2)
library(fitdistrplus)
## Loading required package: MASS
## Loading required package: survival
library(FAdist)
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")
source("cost_functions.R")
## Loading required package: stats4
## Loading required package: splines
##
## Attaching package: 'VGAM'
## The following objects are masked from 'package:FAdist':
##
## dgev, dgumbel, pgev, pgumbel, qgev, qgumbel, rgev, rgumbel
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)))
}
# 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))
plot(x=sort(times), y=yi, log="x",xlim=range(times), type='b',
main=paste("Failure Rate Distribution"),
ylab="Failure Cumulative Distribution [%]",
xlab="Rate [Days]",yaxt="n")
seqn = 1
if (length(yi)>11) seqn =2
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.9,las=1)
legend(x="topleft",legend = c("Empirical","Weibull 3-par fit"), col=c('black','blue'),pch=c(1,2))
t0=find.t0(times)
t1=NULL
#if (t0>0 & t0<min(times))
t1 = times #- 17167 #("2017-01-01")
#else t1=times
#print(length(t1))
q=fitdist(data = t1 ,distr = "weibull")
fit=q$estimate
sd=q$sd
names(t0)="t0"
curve(expr = log(log((1-pweibull3(x,fit[1],fit[2],t0))^-1)), xlim = c(max(times),min(times)), add = T, col="blue", lwd=2)
curve(expr = log(log((1-pweibull3(x,fit[1]+sd[1],fit[2]+sd[2],t0))^-1)), xlim = c(max(times),min(times)), add = T, col="red", lwd=0.5, lty='dashed')
curve(expr = log(log((1-pweibull3(x,fit[1]-sd[1],fit[2]-sd[2],t0))^-1)), xlim = c(max(times),min(times)), add = T, col="red", lwd=0.5, lty='dashed')
hz=h3(times,fit[1],fit[2],t0)
d=data.frame(time=times,hz=hz,shape=fit[1],scale=fit[2],t0,shape.sd=q$sd[1],scale.sd=q$sd[2],ll=q$loglik)
print(fit)
print(ceiling(t0))
#print(ggplot()+xlab("Rate [days]")+ylab("Hazard Rate")+
# ylim(range(d$hz))+
# ggtitle(label = "Hazard Rate Function")+
# geom_point(data = d, mapping = aes(x = times,y = hz)) +
# geom_smooth(data=d,mapping = aes(x = times,y = hz),formula = y~pweibull3(x,fit[1],fit[2],t0), method="loess"))
return(d)
}
au0=build.au0()
l=build.WBS.code(au0 = au0)
code=l$code
au0=l$au0
df=NULL
#c=code$code[20]
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,]
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)
plot(x=m$start_date[2:n], y=m$dt[2:n], xlim=range,
xast="n", type='b',main=paste("Failure rate vs Time for",c),
ylab="Rate [days]",xlab="Calender Time")
axis(side = 1,at = m$start_date[2:n], tick = T,
labels = as.character(m$start_date[2:n],format("%d-%m")), las=3, cex = 0.7)
m$start_date= (as.numeric(m$start_date)) #-tz
fr = data.frame(time = (m$start_date) , rate = m$dt)
fr=fr[2:n,]
d = data.frame(ix=ix,code=as.character(c),hz.plot(fr$rate))
df= rbind(df,d)
}
## [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
## t0
## 7
## [1] "code ix: 8"
## [1] "code: auto1חשמל"
## shape scale
## 1.015244 45.898537
## t0
## 8
## [1] "code ix: 1"
## [1] "code: auto1מכונות"
## shape scale
## 1.018633 12.265879
## t0
## 1
## [1] "code ix: 16"
## [1] "code: auto1מסגרות"
## [1] "code ix: 7"
## [1] "code: auto2בקרות"
## shape scale
## 0.5363409 58.5198257
## t0
## 1
## [1] "code ix: 15"
## [1] "code: auto2חשמל"
## shape scale
## 0.8662477 52.5471828
## t0
## 14
## [1] "code ix: 5"
## [1] "code: auto2מכונות"
## shape scale
## 1.096798 23.777430
## t0
## 1
## [1] "code ix: 10"
## [1] "code: auto2מסגרות"
## shape scale
## 1.978005 63.379084
## t0
## 14
## [1] "code ix: 3"
## [1] "code: auto3בקרות"
## shape scale
## 0.5427151 33.2745163
## t0
## 1
## [1] "code ix: 6"
## [1] "code: auto3חשמל"
## shape scale
## 0.950746 43.239959
## t0
## 4
## [1] "code ix: 45"
## [1] "code: auto3כללי"
## [1] "code ix: 4"
## [1] "code: auto3מכונות"
## shape scale
## 1.070965 18.146961
## t0
## 1
## [1] "code ix: 18"
## [1] "code: auto3מסגרות"
## shape scale
## 0.4367336 60.5920570
## t0
## 1
## [1] "code ix: 9"
## [1] "code: auto4בקרות"
## shape scale
## 0.7747975 46.0443660
## t0
## 1
## [1] "code ix: 11"
## [1] "code: auto4חשמל"
## shape scale
## 1.083907 45.161233
## t0
## 1
## [1] "code ix: 13"
## [1] "code: auto4מכונות"
## shape scale
## 0.8634671 32.6702019
## t0
## 1
## [1] "code ix: 12"
## [1] "code: auto4מסגרות"
## shape scale
## 11.29599 121.39137
## t0
## -79239
## [1] "code ix: 35"
## [1] "code: bullmer0בקרות"
## shape scale
## 0.3443013 245.3746863
## t0
## 5
## [1] "code ix: 34"
## [1] "code: bullmer0חשמל"
## shape scale
## 0.3969858 129.5624038
## t0
## 4
## [1] "code ix: 36"
## [1] "code: bullmer0מכונות"
## shape scale
## 0.3757584 333.4614369
## t0
## 8
## [1] "code ix: 44"
## [1] "code: c-scan1בקרות"
## [1] "code ix: 47"
## [1] "code: c-scan1חשמל"
## shape scale
## 2.195374 21.867571
## t0
## 7
## [1] "code ix: 48"
## [1] "code: c-scan1כללי"
## [1] "code ix: 43"
## [1] "code: c-scan1מכונות"
## shape scale
## 1.139844 54.802070
## t0
## 6
## [1] "code ix: 42"
## [1] "code: c-scan1מסגרות"
## shape scale
## 0.7294383 80.3068103
## t0
## 10
## [1] "code ix: 21"
## [1] "code: eagle0בקרות"
## shape scale
## 0.8799806 47.4798581
## t0
## 5
## [1] "code ix: 19"
## [1] "code: eagle0חשמל"
## shape scale
## 4.219727 118.237821
## t0
## 79
## [1] "code ix: 20"
## [1] "code: eagle0מכונות"
## shape scale
## 0.9071844 20.9649376
## t0
## 1
## [1] "code ix: 31"
## [1] "code: eastman0בקרות"
## shape scale
## 1.064609 15.367798
## t0
## 2
## [1] "code ix: 33"
## [1] "code: eastman0חשמל"
## [1] "code ix: 32"
## [1] "code: eastman0מכונות"
## shape scale
## 0.6886464 34.6784397
## t0
## 1
## [1] "code ix: 41"
## [1] "code: matec0בקרות"
## shape scale
## 1.32875 36.47529
## t0
## 6
## [1] "code ix: 38"
## [1] "code: matec0חשמל"
## shape scale
## 0.3280408 210.0942056
## t0
## 1
## [1] "code ix: 37"
## [1] "code: matec0מכונות"
## shape scale
## 0.8299869 55.3149005
## t0
## 1
## [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
## t0
## 11
## [1] "code ix: 24"
## [1] "code: matrix1חשמל"
## shape scale
## 0.9679526 56.7467799
## t0
## 5
## [1] "code ix: 25"
## [1] "code: matrix1כללי"
## shape scale
## 1.334752 70.536945
## t0
## 5
## [1] "code ix: 22"
## [1] "code: matrix1מכונות"
## shape scale
## 1.03872 13.80723
## t0
## 1
## [1] "code ix: 23"
## [1] "code: matrix1מסגרות"
## shape scale
## 0.8363598 127.3772728
## t0
## 5
## [1] "code ix: 46"
## [1] "code: matrix2בקרות"
## [1] "code ix: 30"
## [1] "code: matrix2חשמל"
## shape scale
## 0.6722939 30.2753450
## t0
## 1
## [1] "code ix: 27"
## [1] "code: matrix2כללי"
## shape scale
## 0.9608599 55.4597957
## t0
## 4
## [1] "code ix: 29"
## [1] "code: matrix2מכונות"
## shape scale
## 0.7587725 15.6259006
## t0
## 1
## [1] "code ix: 28"
## [1] "code: matrix2מסגרות"
## shape scale
## 0.9373358 17.8077947
## t0
## 2
f = file("fr vs time v2.csv","wt")
write.csv(df,f)
close(f)