Section 4.3_MLE_Page 9

### Maximum Likelihood Estimates for the Weibull Distribution using formula
set.seed(1234)

### Generate random values for Weibull Distribution
alpha<- 1.5   # shape
beta<- 1    # scale	
n<- 50     # population
y<- rweibull(n, alpha, beta)
#Sample value y<- c(0.440675384148230,0.365920884094773,1.05815503189297,0.689501584874580,0.290213987509293,0.654138744293414,0.770855902101362,0.292866735507712,0.524718035196492,1.71255576267005,0.780274600426884,2.50073174330735,0.872560946876537,0.211720893670516,2.12992184490308,0.147013582195010,1.43148046875305,0.996330336741661,1.35860524943696,1.76643976410942,0.888580246495540,1.17994271735502,1.73940457639051,0.784847121487458,1.26995968773832,1.03369233214526,0.697818484701977,0.712451567986963,0.907473498515013,0.968561564054612,0.977094167470982,1.17296838945872,1.43292758721935,1.61969142483500,2.58440750125892,0.431628543052501,0.775773323716324,0.473263369693599,0.420820859725220,1.19920860404057,1.39173074658176,1.96121412482026,0.900682067314606,0.500055737452972,0.0807362222672826,0.707392068045203,2.13882254087087,0.663968017136124,1.17814303203075,0.720935190655982)

### Log-likelihood function using formula
vec<-y
weibull_loglik=function(parm){
n<-length(vec)
gamma<-parm[1]
lambda<-parm[2]
loglik<-n*gamma*log(lambda)+n*log(gamma)+(gamma-1)*sum(log(vec))-sum((vec*lambda)^gamma) # Tutorial Equ. 4.6 page 9
return(-loglik)}
weibull<-nlm(weibull_loglik,parm<-c(1,1), hessian=TRUE) # Tutorial Equ. 4.7, 4.8 for differentiation page 9

### Generate empirical distribution using failure counts and time units
l<-c(seq(0, 3, by = 0.061))  # Random failure counts
p<-c(seq(0.02, 1, by = 0.020)) # Random time units within the desired range

### Cumulative density function plot
plot(vec,p,xlab = "Time (t)", 
      ylab = "Cumulative density function", main = "Fitted Weibull density using MLE")
par(new=TRUE)

plot(function(l) pweibull(l, shape = weibull$estimate[1], 
      scale = 1/weibull$estimate[2]), 0, 3, axes=False)

### Probability density function plot
hist(vec,xlab = "Time (t)", 
      ylab = "Probability density function", main = "Fitted Weibull density using MLE")
par(new=TRUE)
plot(function(vec) dweibull(vec, shape = weibull$estimate[1], 
      scale = 1/weibull$estimate[2]), 0, 3, axes=False)