#!/usr/local/bin/pythonw
# coding: iso-8859-1
# Copyright © 2009 by Amos Newcombe
# 
# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License along with this program.  If not, see <http://www.gnu.org/licenses/>.

'''the Beta function: complete and incomplete
'''

from math import exp, log, fabs
from Gamma import LnGamma

#	Beta function

def Beta(a, b):
	return exp(LnGamma(a) + LnGamma(b) - LnGamma(a+b))

#	Incomplete Beta function

def IBeta(a, b, x):
	if x  <= 0.: return 0.0
	if 1. <= x : return 1.0
	a, b = float(a), float(b)
	#if x == 0.0 or x == 1.0: return x
	xFactor = exp(LnGamma(a+b) - LnGamma(a) - LnGamma(b) + a*log(x) + b*log(1.0-x))
	if x < (a+1.0)/(a+b+2.0):
		xResult = xFactor * IBetaCF(a, b, x) / a
	else:
		xResult = 1.0 - xFactor * IBetaCF(b, a, 1.0-x) / b
	return xResult

def IBetaCF(a, b, x):
	iMax = 101
	epsilon = 3.0e-7
	xTiny = 1.0e-30
	xAB = a+b
	xAP = a+1.0
	xAM = a-1.0
	c = 1.0
	d = 1.0 - x*xAB/xAP
	if fabs(d) < xTiny: d = xTiny
	d = 1.0/d
	h=d
	i = 1
	while i < iMax:
		i2 = 2*i
		aa = i * (b-i) * x / ((xAM+i2) * (a+i2))
		d = 1.0 + aa*d
		if fabs(d) < xTiny: d = xTiny
		c = 1.0 + aa/c
		if fabs(c) < xTiny: c = xTiny
		d = 1.0/d
		h = h * d * c
		aa = -(a+i) * (xAB+i) * x / ((xAP+i2) * (a+i2))
		d = 1.0 + aa*d
		if fabs(d) < xTiny: d = xTiny
		c = 1.0 + aa/c
		if fabs(c) < xTiny: c = xTiny
		d = 1.0/d
		xDC = d * c
		h = h * xDC
		if fabs(xDC - 1.0) < epsilon:
			return h
		i = i+1
	raise ValueError, 'IBetaCF('+`a`+', '+`b`+', '+`x`+') fails'

def InverseIBeta(a, b, y, xEpsilon=1e-7):
	if y == 0.: return 0.
	if y == 1.: return 1.
	a, b = float(a), float(b)
	x = a / (a+b)
	iMax = 101
	i = 1
	while i < iMax:
		d = y - IBeta(a, b, x)
		if fabs(d) < xEpsilon: return x
		m = pow(x, a-1.0) * pow(1.0-x, b-1.) / Beta(a, b)
		#print x, d, m
		x = x + d / m
		if 1. < x: x = 1. - xEpsilon
		if x < 0.: x = xEpsilon
	raise ValueError, 'InverseIBeta('+`a`+', '+`b`+', '+`x`+') fails'

if __name__ == '__main__':
	pass