# loglikelihood function

log.like<-function(x,a,b,g) {
length(x)*log(b)-length(x)*log(gamma(1+a))+sum(log(gamma(x+a)))+log(g)*sum(x-1))
}

#my loglikelihood
log.like<-function(x,a,b,g) {
length(x)*log(b)-length(x)*log(gamma(1+a))+sum(log(gamma(x+a)))+log(g)*(sum(x)-length(x))
}

# prepare the data

a<-read.table(file="table1.8.dat")
s<-1
y<-NULL
for(i in 1:10) {for(j in 1:10) {y[s]<-a[i,j]; s<-s+1}}
s<-1
x<-NULL
for(i in 1:100) {for(j in 1:y[i]) {x[s]<-i; s<-s+1}}
str(x)

# initial values

dx<-0.01

a<-0.1
b<-2.5
g<-0.01

etemp<-NULL

for (i in 1:50) {

f0<-log.like(x,a,b,g)

f1<-log.like(x,a+dx,b,g)
f2<-log.like(x,a,b+dx,g)
f3<-log.like(x,a,b,g+dx)

b1<-(f1-f0)/dx
b2<-(f2-f0)/dx
b3<-(f3-f0)/dx

f11<-log.like(x,a+2*dx,b,g)
f12<-log.like(x,a+dx,b+dx,g)
f13<-log.like(x,a+dx,b,g+dx)
f22<-log.like(x,a,b+2*dx,g)
f23<-log.like(x,a,b+dx,g+dx)
f33<-log.like(x,a,b,g+2*dx)

g11<-(f11-2*f1+f0)/dx^2
g22<-(f22-2*f2+f0)/dx^2
g33<-(f33-2*f3+f0)/dx^2
g12<-(f12-f1-f2+f0)/dx^2
g13<-(f13-f1-f3+f0)/dx^2
g23<-(f23-f2-f3+f0)/dx^2

est<-rbind(a,b,g)

bb<-cbind(b1,b2,b3)
gg<-cbind(c(g11,g12,g13),c(g12,g22,g23),c(g13,g23,g33))

etemp<-est-solve(gg)%*%t(bb)

a<-etemp[1]
b<-etemp[2]
g<-etemp[3]

cat("Iteration ",i," ",round(a,3)," ",round(b,3)," ",round(g,3)," ","\n")

if(sum(abs(bb))<0.00001) break

}