set.seed(432)
n<-6
X<-rweibull(n,2)
X<-round(X,3)

w.loglike<-function(av,xv){

	m1<-t(outer(av-1,log(xv),"*"))-outer(xv,av,"^")+t(outer(log(av),rep(1,length(xv)),"*"))
	return(apply(m1,2,sum))

}

w.diffloglike<-function(av,xv){

	m1<-t(outer(rep(1,length(av)),log(xv),"*"))-outer(xv,av,"^")*outer(log(xv),rep(1,length(av)),"^")+
	    t(outer(1/av,rep(1,length(xv)),"*"))
	return(apply(m1,2,sum))

}

w.ddiffloglike<-function(av,xv){

	m1<--outer(xv,av,"^")*outer((log(xv))^2,rep(1,length(av)),"^")-
	    t(outer(1/av^2,rep(1,length(xv)),"*"))
	return(apply(m1,2,sum))

}


avec<-seq(0,4,by=0.001)

fvec<-w.loglike(avec,X)
dfvec<-w.diffloglike(avec,X)
ddfvec<-w.ddiffloglike(avec,X)


pdf("Weibull-loglike.pdf",paper="USr",width=11,height=9)
plot(avec,fvec,type="l",ylab=expression(l(theta)),xlab=expression(alpha))
dev.off()
pdf("Weibull-dloglike.pdf",paper="USr",width=11,height=9)
plot(avec,dfvec,type="l",xlim=range(0,4),ylim=range(-10,10),
     ylab=expression(dot(l)(theta)),xlab=expression(alpha));abline(h=0,col="red")
dev.off()
pdf("Weibull-ddloglike.pdf",paper="USr",width=11,height=9)
plot(avec,ddfvec,type="l",ylim=range(-20,-2),
     ylab=expression(ddotl(theta)),xlab=expression(alpha));abline(h=0,col="red")
dev.off()

########################################################################

x0<-0.25
d0<-w.diffloglike(x0,X)
dd0<-w.ddiffloglike(x0,X)
a0<-d0-dd0*x0
plot(avec,dfvec,type="l",ylim=range(-30,30),ylab="dg(x)");abline(h=0,col="red")
points(x0,d0,pch=4,lwd=2)
abline(a0,dd0,lty=2)

x1<-x0-d0/dd0
points(x1,0,pch=19,cex=2)
d1<-w.diffloglike(x1,X)
dd1<-w.ddiffloglike(x1,X)
a1<-d1-dd1*x1
points(x1,d1,pch=4,lwd=2)
abline(a1,dd1,lty=2)

x2<-x1-d1/dd1
d2<-w.diffloglike(x2,X)
dd2<-w.ddiffloglike(x2,X)
a2<-d2-dd2*x2
points(x2,d2,pch=4,lwd=2)
abline(a2,dd2,lty=2)

x3<-x2-d2/dd2
d3<-w.diffloglike(x3,X)
dd3<-w.ddiffloglike(x3,X)
a3<-d3-dd3*x3
points(x3,d3,pch=4,lwd=2)
abline(a3,dd3,lty=2)

x4<-x3-d3/dd3
d4<-w.diffloglike(x4,X)
dd4<-w.ddiffloglike(x4,X)
a4<-d4-dd4*x4
points(x4,d4,pch=4,lwd=2)
abline(a4,dd4,lty=2)


###########################################################

xcur<-0.25
dcur<-w.diffloglike(xcur,X)
ddcur<-w.ddiffloglike(xcur,X)
acur<-dcur-ddcur*xcur
plot(avec,dfvec,type="l",xlab=expression(alpha),ylim=range(-30,30),ylab="dg(x)");abline(h=0,col="red")
points(xcur,dcur,pch=4,lwd=2)
abline(acur,ddcur,lty=2)
fname<-paste("Run1-Step-",0,".pdf",sep="")
dev.copy2pdf(file=fname,height=9,width=11)
for(iter in 1:6){	
	fname<-paste("Run1-Step-",iter,".pdf",sep="")
	xcur<-xcur-dcur/ddcur
	points(xcur,0,pch=19,cex=1.2)
	dcur<-w.diffloglike(xcur,X)
	ddcur<-w.ddiffloglike(xcur,X)
	acur<-dcur-ddcur*xcur
	points(xcur,dcur,pch=4,lwd=2)
	abline(acur,ddcur,lty=2)
	abline(v=xcur,lty=2)
	text(2.0,-15-2*iter,paste(iter,": ",round(xcur,4)),cex=1.5)
	dev.copy2pdf(file=fname,height=9,width=11)
	print(c(iter,xcur,dcur))
}

###########################################################

xcur<-3.25
dcur<-w.diffloglike(xcur,X)
ddcur<-w.ddiffloglike(xcur,X)
acur<-dcur-ddcur*xcur
plot(avec,dfvec,type="l",xlab=expression(alpha),ylim=range(-30,30),ylab="dg(x)");abline(h=0,col="red")
points(xcur,dcur,pch=4,lwd=2)
abline(acur,ddcur,lty=2)
fname<-paste("Run2-Step-",0,".pdf",sep="")
dev.copy2pdf(file=fname,height=9,width=11)
for(iter in 1:6){	
	fname<-paste("Run2-Step-",iter,".pdf",sep="")
	xcur<-xcur-dcur/ddcur
	points(xcur,0,pch=19,cex=1.2)
	dcur<-w.diffloglike(xcur,X)
	ddcur<-w.ddiffloglike(xcur,X)
	acur<-dcur-ddcur*xcur
	points(xcur,dcur,pch=4,lwd=2)
	abline(acur,ddcur,lty=2)
	abline(v=xcur,lty=2)
	text(2.0,-15-2*iter,paste(iter,": ",round(xcur,4)),cex=1.5)
	dev.copy2pdf(file=fname,height=9,width=11)
	print(c(iter,xcur,dcur))
}

#############################################################################################
#Bisection

aval<-1
bval<-3
xval<-(aval+bval)/2
avec<-bvec<-xvec<-ivec<-numeric()

avec<-aval
bvec<-bval
ivec<-0
xvec<-xval


for(iter in 1:20){

	d1<-w.diffloglike(aval,X)
	d2<-w.diffloglike(xval,X)
	if(d1*d2 <= 0){
		bval<-xval
	}else{
		aval<-xval
	}
	xval<-(aval+bval)/2

	avec<-c(avec,aval)
	bvec<-c(bvec,bval)
	ivec<-c(ivec,iter)
	xvec<-c(xvec,xval)
	print(c(iter,round(c(aval,bval,xval),5)))
}
data.frame(iter=ivec,a=avec,b=bvec,x=xvec)

library(survival)

survreg(Surv(X)~-1,dist="weibull")

mean(X)