lifetimes<-c(1051,1337,1389,1921,1942,2322,3629,4006,4012,4063, 4921,5445,5620,5817,5905,5956,6068,6121,6473,7501,7886,8108,8546,8666,8831,9106,9711,9806,10205,10396, 10861,11026,11214,11362,11604,11608,11745,11762,11895,12044,13520,13670,14110,14496,15395,16179,17092,17568,17568)




x=seq(0,10, length=100)
par(mfrow=c(1,2))
plot(x, dweibull(x, shape=2, scale=1), type="l", ylab="Density")
lines(x, dweibull(x, shape=2, scale=2), lty=2)
title(main="Weibull Distribution, shape=2")
leg.names<-c("scale=1","scale=2")
legend(locator(1),leg.names,lty=1:2,bty="n")



loglikelihood <- function(data, theta, lambda=2)
 { logl<-sum((lambda-1)*log(data)+log(lambda)-lambda*log(theta)-(data/theta)^{lambda})
   return(logl)
 }


theta1 <- seq(7000, 15000, by=100)
 loglik=rep(NA, length(theta1))
 for (i in 1:length(theta1))
 {
 loglik[i]=loglikelihood(lifetimes, theta1[i])
 }


F <- function(data, theta, lambda=2)
 {fx<--(lambda*length(data)/theta)+
 lambda*(sum(data^{lambda}))/(theta^{3})
 return(fx)
 }

devF<-function(data, theta)
 {dfx<-2*length(data)/(theta^2)-6*(sum(data^2))/(theta^{4})
 return(dfx)
 }

theta0<-8805.694
test<-3
while (test!=10) {
theta.new<-theta0-F(lifetimes,theta0,lambda=2)/devF(lifetimes,theta0)

if ((theta.new-theta0)<=10^(-6)) {
                 test<-10} else {
                     print(theta.new-theta0)
                     theta0<-theta.new
                      }
}