set.seed(432)
n<-20
X<-rweibull(n,2)

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)


plot(avec,fvec,type="l")
plot(avec,dfvec,type="l",xlim=range(1,3),ylim=range(-10,10),ylab="dg(x)");abline(h=0,col="red")
plot(avec,ddfvec,type="l",ylab="ddg(x)",ylim=range(-20,-2));abline(h=0,col="red")

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

x1<-x0-d0/dd0
d1<-w.diffloglike(x1,X)
dd1<-w.ddiffloglike(x1,X)
a1<-d1-dd1*x1
points(x1,d1,pch=4,lwd=5)
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=5)
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=5)
abline(a3,dd3,lty=2)









avec[which.max(fvec)]