# y_scale_coeff is the
# coefficient to reduce y axis values
# only to display
y_scale_coeff = 1/1000

min_time = 0
max_time = 10000
max_time = max_time * y_scale_coeff
xlim_bottom = c(0,200)

nx = 20
ny = 40

ylim = c(min_time,max_time)

# x
# Cum. Hazard % values for faiure hours
# diesel engine fan NELSON data
x = c(1.4,2.9,4.4,5.9,7.7,9.6,11.5,13.6,15.8,18.7,22.7,35.2)

# y
# faiure hours
# diesel engine fan NELSON data
y = c(450,1150,1150,1600,2070,2070,2080,3100,3450,4600,6100,8750)

y = y * y_scale_coeff

# ticks at bottom x axis
ticks_x = c(0,50,100,150,200)

# ticks at y axis
ticks_y = seq(0, max_time, 1)

# labels at top x axis
x_top_labels = 100 * (1 - 1 / (exp(1) ^ (ticks_x / 100)))
x_top_labels = round(x_top_labels,digits = 1)

# plot
plot(x,y,bty="n",xlim=xlim_bottom,ylim=ylim,xaxs="i",yaxs="i",col="red",pch=16,axes=FALSE,xlab="",ylab="")

#x axis at bottom
axis(1,ticks_x,col="red",col.ticks="red",col.axis="red")
mtext("Cum. Hazard %",1,line=2.5,at=100,col="red")

#x axis at top
axis(3,ticks_x,labels=x_top_labels,col="blue",col.ticks="blue",col.axis="blue")
mtext("CDF %",3,line=2.5,at=100,col="blue")

# y axis on left
axis(2,ticks_y,col="black",col.ticks="black",col.axis="black")
mtext(paste("Failure Hours ( X ",y_scale_coeff,")"),2,line=2.5,at=max_time/2,col="black")

# add grid to the plot
grid(nx,ny,col="lightgray",lty="solid",lwd=par("lwd"),equilogs=TRUE)

# add points again
# (to layer dots over grid lines)
points(x,y,col="red")

# draw line to determine mttf
# vertical line from top to Cum Haz. % = 100
segments(100,0,100,max_time,col="green")

# get best linear fit
fit = lm(y~x)

# draw best fit line
abline(fit)

summary(fit)

newdata = data.frame(x=100)
predict(fit,newdata,interval="predict")
## R prediction 24538.7
## Wayne Nelson's prediction 29000
## Best Numeric Prediction 28703.3