############################
#
#  Disequilibrium coefficient
#
###########################


set.seed(1234)			## seed
m<-10000000			## Monte Carlo sample size

a1<-0.5				##
a2<-0.5				## Prior parameters
a3<-0.5				##

n1<-4				##
n2<-18				## Observed frequencies
n3<-94				##

p1<-rgamma(m,a1+n1,2)
p2<-rgamma(m,a2+n2,2)
p3<-rgamma(m,a3+n3,2)
theta<-4*p1*p3/(p2^2)
lambda<- (theta^0.5 - 1)/(theta^0.5 + 1)
lambda<-c(-1,lambda,1)
freq<-hist(lambda,breaks=200,main="Posterior of lambda",probability=T,ylab="Density")
count_pos<-sum(freq$count[1:100])
position_pos<-200
while (position_pos>=101){
   if (freq$count[position_pos]<=freq$count[100])count_pos<-count_pos+freq$count[position_pos]
   position_pos<-position_pos-1
}
count_neg<-sum(freq$count[101:200])
position_neg<-1
while (position_neg<=100){
   if (freq$count[position_neg]<=freq$count[101])count_neg<-count_neg+freq$count[position_neg]
   position_neg<-position_neg+1
}
evalue_pos<-(count_pos-2)/m
evalue_neg<-(count_neg-2)/m

#### HPD - Credible Interval ####
alpha<-0.05				## (1-alpha) is the credibility
n.HPD<-m*alpha
LI<-1
LS<-200
count.HPD<-freq$count[LI]+freq$count[LS]
while (count.HPD < n.HPD){
   if (freq$count[LI]<=freq$count[LS])
      {LI<-LI+1
       count.HPD<-count.HPD+freq$count[LI]
      }
   if (freq$count[LI]> freq$count[LS])
      {LS<-LS-1
       count.HPD<-count.HPD+freq$count[LS]
      }
}

mean(lambda)		## bayes estimate of lambda

(LI-100)/100		## Lower Limit of Credible Interval
(LS-100)/100		## Upper Limit of Credible Interval

evalue_pos		## e-value, if lambda is positive
evalue_neg		## e-value, if lambda is negative