#modifications by John Christie
#03/10/2011
# The variable "slices" in each function is there so that one can adjust the quality of the measure of the area under the curve.  It's one of the main reasons that this code will give different results from online calculators.  It can also be lowered in order to overcome any possible performance issues (unwise in the cross_2… functions).    Everything should be fast enough as is.

Bf<-function(sd, obtained, uniform = FALSE, lower=0, upper=1, meanoftheory=0, sdtheory=1, tails=2)
{
#Version 2.0.1
#modifications by John Christie
# modification means that this does not exactly replicate Dienes.  What it does
# is do the right thing instead.  :)  The current version is more accurate.
#03/10/2011
#Authors Danny Kaye & Thom Baguley
#Version 1.0
#19/10/2009
# test data can be found starting at p100
# notes on how to use these functions can be found at
# http://danny-kaye.co.uk/Docs/Dienes_notes.pdf
	#raised from 2000 for better accuracy - the speed of the new code allows it
	slices <- 20000
	if(uniform){
		range <- upper - lower
		dist_theta <- 1 / range
		# incr <- range / slices
		# theta <- seq(lower + incr, by = incr, length.out = slices+1)
		# height <- dist_theta * dnorm(obtained, theta, sd)
		# area <- sum(height * incr)
		# the commented code above replicates the original result but at the
		# limit (slices <- 5e6) it's actually equivalent to the following…
		area <- dist_theta * diff(pnorm(c(lower, upper), obtained, sd)) 
	}else{
		# the code below again doesnt' replicate the original code.
		# incrementing in a scalar loop was causing an accumulation of tiny fp errors
		# the lower end was incremented prior to starting (that's fixed above too)
		zlim <- 5
		incr <- sdtheory / (slices/(zlim*2))
		newLower <- meanoftheory - zlim * sdtheory
		theta <- seq(newLower, by = incr, length.out = slices+1)
		dist_theta <- dnorm(theta, meanoftheory, sdtheory)
		if (tails == 1){
			dist_theta <- dist_theta[theta > 0]	* 2
			theta <- theta[theta > 0]	
			}
		height <- dist_theta * dnorm(obtained, theta, sd)
		area <- sum(height * incr)
	}
	LikelihoodTheory <- area
	Likelihoodnull <- dnorm(obtained, 0, sd)
	BayesFactor <- LikelihoodTheory / Likelihoodnull
	return( list("LikelihoodTheory" = LikelihoodTheory, "Likelihoodnull" = Likelihoodnull, "BayesFactor" = BayesFactor) )
}


contrast<-function(SEd, meand, nu, theta1, theta2)
{
#modifications by John Christie
#01/10/2011
#Authors Danny Kaye & Thom Baguley
#19/10/2009
#The test data is found in Box 3.7 pp 74-5
#SEd - contrast standard error
#meand - contrast mean
#nu - contrast  degrees of freedom
#
	slices <- 10000
	zlim <- 5
	inc <- SEd / (slices/(zlim*2))

	theta <- seq(meand - zlim * SEd, by = inc, length.out = slices+1)
	likelihood <- (1 + (meand - theta)^2 / (nu * SEd^2))^(-(nu + 1) / 2)
	likelihoodmax <- max(likelihood)

	likelihood <- likelihood / likelihoodmax

	which32 <- which(likelihood >= 1/32)
	which8 <- which(likelihood >= 1/8)
	begin32 <- theta[which32[1]]
	begin8 <- theta[which8[1]]
	end8 <- theta[tail(which8,1)+1]
	end32 <- theta[tail(which32,1)+1]

	B1 <- round((theta1 - theta[1]) / inc + 1, 0)
	B2 <- round((theta2 - theta[1]) / inc + 1, 0)
	likelihoodratio <- likelihood[B1] / likelihood[B2]
	return( list("begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32, "likelihoodratio"=likelihoodratio) )
}

meanUV<-function(SEd, meand, n, theta1, theta2)
{
#modifications by John Christie
#01/10/2011
#Authors Danny Kaye & Thom Baguley
#Version 1.0
#19/10/2009
#The test data is from page 138
#SEd - sample standard error
#meand - sample mean
#n - number of subjects
#theta1 - population mean assumed by first hypothesis
#theta2 - population mean assumed by second hypothesis
#
	Vard <- n * SEd^2
	SSd <- Vard * (n - 1)

	slices <- 10000
	zlim <- 5
	inc <- SEd / (slices/(zlim*2))

	theta <- seq(meand - zlim * SEd, by = inc, length.out = slices+1)
	likelihood <- (SSd + n * (meand - theta)^2)^(-(n - 2) / 2)
	likelihoodmax <- max(likelihood)

	likelihood <-likelihood / likelihoodmax

	#unlike the original I find a chunk that's the middle of the distribution
	which32 <- which(likelihood >= 1/32)
	which8 <- which(likelihood >= 1/8)
	begin32 <- theta[which32[1]]
	begin8 <- theta[which8[1]]
	end8 <- theta[ tail(which8, 1) + 1 ]
	end32 <- theta[ tail(which32, 1) + 1 ]

	B1 <- round((theta1 - theta[1]) / inc + 1, 0)
	B2 <- round((theta2 - theta[1]) / inc + 1, 0)
	likelihoodratio <- likelihood[B1] / likelihood[B2]
	return( list("begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32, "likelihoodratio"=likelihoodratio) )
}


cross2x2_odds <- function(a, b, c, d, psi1, psi2)
{
#modifications by John Christie
#01/10/2011
#Authors Danny Kaye & Thom Baguley
#Version 1.0
#19/10/2009
#The test data is from figure 5.4 p 135
#
	m <- a + b
	n <- c + d
	k <- a + c
	l <- b + d

	slices <- 10000
	psi <- (1:slices) / sqrt(slices)

	starting <- 0
	if((a - d) > 0) starting <- a - d
	ending <- m
	if((a + c) < m) ending <- a + c

	A <- starting:ending
	likelihood <- sapply( 1:slices, function(B) {
			sum( choose(m, A) * choose(n, (a + c - A)) * psi[B]^(A - a))})

	likelihood <- likelihood^-1
	likelihoodmax <- max(likelihood)
	psimax <- psi[which.max(likelihood)]

	likelihood <- likelihood / likelihoodmax

	which32 <- which(likelihood >= 1/32)
	which8 <- which(likelihood >= 1/8)
	begin32 <- psi[which32[1]]
	begin8 <- psi[which8[1]]
	end8 <- psi[ tail(which8, 1) + 1 ]
	end32 <- psi[ tail(which32, 1) + 1 ]
	
	psi1 <- psi1 * 100
	psi2 <- psi2 * 100
	likelihoodratio <- likelihood[psi1] / likelihood[psi2]
	return( list("psimax"=psimax, "begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32, "likelihoodratio"=likelihoodratio) )
}

cross2x22<-function(a, b, c, d, gamma1, gamma2)
{
#modifications by John Christie
#01/10/2011
#Authors Danny Kaye & Thom Baguley
#Version 1.0
#19/10/2009
#The test data is from figure 5.4 p 135
#
	m <- a + b
	n <- c + d
	k <- a + c
	l <- b + d

	slices <- 10000
	gammamax <- 0
	gamma <- (1:slices) / sqrt(slices)

	A <- b:(m + n - d)
	likelihood <- sapply( 1:slices, function(B){
		sum( choose(A - 1, b - 1) * choose(m + n - A - 1, d - 1) * gamma[B]^(A - m) ) } )

	likelihood <- likelihood^-1
	likelihoodmax <- max(likelihood)
	gammamax <- gamma[which.max(likelihood)]

	likelihood <- likelihood / likelihoodmax

	which32 <- which(likelihood >= 1/32)
	which8 <- which(likelihood >= 1/8)
	begin32 <- gamma[which32[1]]
	begin8 <- gamma[which8[1]]
	end8 <- gamma[ tail(which8, 1) + 1 ]
	end32 <- gamma[ tail(which32, 1) + 1 ]
	
	gamma1 <- gamma1 * 10
	gamma2 <- gamma2 * 10
	likelihoodratio <- likelihood[gamma1] / likelihood[gamma2]
	return( list("gammamax"=gammamax, "begin8"=begin8, "end8"=end8, "begin32"=begin32, "end32"=end32, "likelihoodratio"=likelihoodratio) )
}

proportion <- function(suc, fail)
{
#modifications by John Christie
#01/10/2011
#Authors Danny Kaye & Thom Baguley
#Version 1.0
#19/10/2009
#The test data is from figure 5.3 p131
#not sure this does what he says it does, no theta input is asked for
#in the book code, as implied in the text and no ratio is returned
#note that the function is called with successes and failures not 
#successes and trials
#
	slices <- 10000
	theta <- (1:slices)/slices

	likelihood <- theta^suc*(1-theta)^fail
	likelihoodmax <- max(likelihood)
	thetamax <- theta[which.max(likelihood)]

	likelihood <- likelihood / likelihoodmax

	#unlike the original I find a chunk that's the middle of the distribution
	which32 <- which(likelihood >= 1/32)
	which8 <- which(likelihood >= 1/8)
	begin32 <- theta[which32[1]]
	begin8 <- theta[which8[1]]
	end8 <- theta[ tail(which8, 1) + 1 ]
	end32 <- theta[ tail(which32, 1) + 1 ]

	return( list("thetamax"=thetamax, "begin8"=begin8, "end8"=end8,
	"begin32"=begin32, "end32"=end32) )
}