counts<-rep(0,153)
num<-matrix(counts,nrow=17,ncol=9)
num[1,1]<-1
num[2,1]<-1
num[4,1]<-3
num[5,1]<-1
num[5,2]<-1
num[6,1]<-1
num[6,2]<-2
num[7,1]<-1
num[7,2]<-3
num[7,4]<-1
num[8,1]<-2
num[8,3]<-1
num[9, ]<-c(5,5,2,1,0,0,0,0,0)
num[10,]<-c(17,6,2,0,0,0,0,0,0)
num[11,]<-c(14,9,3,1,1,2,0,0,0)
num[12,]<-c(19,17,7,2,2,0,0,0,0)
num[13,]<-c(9,14,7,3,0,0,2,1,0)
num[14,]<-c(10,5,5,1,1,0,0,0,1)
num[15,]<-c(4,1,2,0,0,0,1,0,0)
num[16,2]<-2
num[16,3]<-1
num[17,1]<-2

prob<-function(n,k,alpha,beta){
	gamma(alpha+k)*gamma(n-k+beta)*gamma(alpha+beta)/(gamma(n+alpha+beta)*gamma(alpha)*gamma(beta))
}

loglik<-function(paras){
	loglike<-0
	alpha<-paras[1]
	beta<-paras[2]
	for(i in 1:17){
		for(j in 1:9){
			if (num[i,j]!=0)
				loglike<-loglike+num[i,j]*log(prob(i,j-1,alpha,beta))
			}
		-loglike
	}
}


#factorial
fac<-function(n){
	fac<-1
	for(i in 1:n){
	fac<-fac*i
}
fac<-ifelse(n==0,1,fac)
}

#combination
comb<-function(n,m){
	com<-fac(n)/(fac(m)*fac(n-m))
	com<-ifelse(n<m,0,com)
	com
}

estimation<-function(estimates) {
	counts<-rep(0, 153)
	est<-matrix(counts,nrow=17)
	for(i in 1:17){ 
		n<-sum (num[i,])
		est[i,1]<-n*prob(i,0, estimates[1], estimates[2])
		for(j in 2:9){
			est[i,j]<-n*prob(i,j-1,estimates[1], estimates[2])*comb(i,j-1)
		}
	}
	est
}

brooks.mle<-nlm(loglik,c(1,1))
brooks.mle
brooks.table<-estimation(c(brooks.mle$estimate[1],brooks.mle$estimate[2]))
round(brooks.table,1)