#####Weibull Distribution

x=seq(0,10, length=100)
plot(x, dweibull(x, shape=2, scale=1), type="l", ylab="Density")
points(x, dweibull(x, shape=2, scale=2), type="l", col=2)

plot(x, dweibull(x, shape=1, scale=1), type="l", ylab="Density")
points(x, dweibull(x, shape=1, scale=2), type="l", col=2)



#####Read the data "Failure Times for pressure levels"

failuretimes <- read.table("Failure times of pressure
vessels.txt",header=T)

failuretimes

summary(failuretimes)

#    lifetimes
# Min.  : 1051
# 1st Qu.: 5620
# Median : 8831
# Mean  : 8806
# 3rd Qu.:11745
# Max.  :17568

attach(failuretimes)

#######Histogram of the Data###########

hist(lifetimes)

###########qq.plot for Weibul Distribution with lambda=2######

library("car")
qq.plot(lifetimes, dist="weibull", shape=2)

qqline(lifetimes, rweibull(49, shape=2))


#############Maximize the likelihood function with
#############Newton Raphson Method

###Define the loglikelihood function

loglikelihood <- function(data, theta, lambda=2)
{
logl <- rep(NA, length(theta))
for (i in 1:length(theta))
{
logl[i] <- sum( (lambda-1)*log(data)+log(lambda)-
            lambda*log(theta[i])-(data/theta[i])^{lambda})
}
return(logl)
}

theta1 <- seq(7000, 15000, by=100)

plot(theta1, loglikelihood(lifetimes, theta1), type="l",
xlab="Values of theta", ylab="log likelihood function")


####Define the score function

get.score <- function(data, theta, lambda=2)
{
score <- rep(NA, length(theta))
for (i in 1:length(theta))
{
score[i] <-
-(lambda*length(data)/theta[i])+2*(sum(data^{2}))/(theta[i]^{3})
}
return(score)
}

##########Define the Information matrix

get.information <- function(data, theta, lambda=2)
{
information <- (lambda^{2}*length(data))/(theta^{2})
return(information)
}

###########Newton Raphson for estimating the
###########parameter theta

ybar <- mean(lifetimes)
ybar
theta <- ybar
it  <- 0  ####iterative count
del <- 1  ####iterative adjustment
while(abs(del) > 0.00001 && (it <- it+1) < 10)
{
del <- get.score(lifetimes, theta)/get.information(lifetimes, theta)
theta <- theta+del
cat(it, theta, get.score(lifetimes, theta), get.information(lifetimes,
theta),"\n")
}

#1 9959.204 -0.0001320064 1.976090e-06
#2 9892.402 -4.517492e-07 2.002869e-06
#3 9892.177 -5.150392e-12 2.002960e-06
#4 9892.177 1.734723e-18 2.002960e-06

sqrt(mean(lifetimes^{2}))

#9892.177  ####exact value

##########Confidence Interval

sderror <- sqrt(1/2.002960e-06)


9892.177-1.96*sderror;    9892.177+1.96*sderror

#[1] 8507.272
#[1] 11277.08