############################
#
#  Yule's coefficient of association
#
###########################


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

a1<-0.5				##
a12<-0.5			##    Prior parameters
a21<-0.5			##
a3<-0.5				##

n1<-28				##
n12<-12				## Observed frequencies
n21<-1				##
n3<-4				##

q1<-rgamma(m,a1+n1,2)
q2<-rgamma(m,a12+n12,2)
q3<-rgamma(m,a21+n21,2)
q4<-rgamma(m,a3+n3,2)
lambda<- ((q1*q4)^.5 - (q2*q3)^.5)/((q1*q4)^.5 + (q2*q3)^.5)
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