n<-500
set.seed(2342)

a<-2
b<-1
X<-rweibull(n,shape=2,scale=1)

log.like.w<-function(abv,Xv){

	tvec<-sum(dweibull(Xv,shape=abv[1],scale=abv[2],log=T))
	return(tvec)

}



log.like.w(c(2,1),X)


avec<-seq(1,3,by=0.02)
bvec<-seq(0.5,1.5,by=0.01)
lmatw<-matrix(0,nrow=length(avec),ncol=length(bvec))
for(i in 1:length(avec)){
	for(j in 1:length(bvec)){
	
		lmatw[i,j]<-log.like.w(c(avec[i],bvec[j]),X)

	}
	
}

image(avec,bvec,lmatw,col=heat.colors(50),xlab=expression(a),ylab=expression(b))
contour(avec,bvec,lmatw,nlevels=50,add=T)
abline(v=a,h=b,col="red")

######################################################################
log.like.w2<-function(abv,Xv){

	av<-exp(abv[1])
	bv<-exp(abv[2])

	tvec<-sum(dweibull(Xv,shape=av,scale=bv,log=T))
	return(tvec)

}
optim(c(0,0),fn=log.like.w2,Xv=X,control=list(fnscale=-1))



logregw.opt<-optim(c(0,0),fn=log.like.w2,Xv=X,
                  method="BFGS",control=list(fnscale=-1),hessian=T)
image(avec,bvec,lmatw,col=heat.colors(50),xlab=expression(a),ylab=expression(b))
contour(avec,bvec,lmatw,nlevels=50,add=T)
abline(v=exp(logregw.opt$par[1]),h=exp(logregw.opt$par[2]),col="red")


Asymp.var.log<-solve(-logregw.opt$hessian)
Asymp.var<-diag(exp(logregw.opt$par)) %*% Asymp.var.log %*% diag(exp(logregw.opt$par))
                         


library(MASS)
Approx.AB<-mvrnorm(2000,mu=exp(logregw.opt$par),Sigma=Asymp.var)
points(Approx.AB,pch=19,cex=0.5)

avec.new<-seq(1.8,2.2,by=0.002)
bvec.new<-seq(0.8,1.2,by=0.002)

lmat2w<-matrix(0,nrow=length(avec.new),ncol=length(bvec.new))
for(i in 1:n){
	lmat2w<-lmat2w+outer(avec.new,bvec.new,function(x,y,Xi){return(dweibull(Xi,shape=x,scale=y,log=T))},Xi=X[i])
}

image(avec.new,bvec.new,lmat2w,col=heat.colors(50),xlab=expression(a),ylab=expression(b))
contour(avec.new,bvec.new,lmat2w,add=T,nlevels=30)
points(Approx.AB,pch=19,cex=0.5)