func genBetaChengBB(aa, bb, a, b float64) float64 {
	alpha := a + b
	beta := math.Sqrt((alpha - 2.0) / (2.0*a*b - alpha))
	gamma := a + 1.0/beta
	var r, s, t, v, w, z float64
	complete := func() bool {
		u1, u2 := rand.Float64(), rand.Float64()
		v = beta * math.Log(u1/(1.0-u1))
		if v <= 709.78 {
			w = a * math.Exp(v)
			if math.IsInf(w, 0) {
				w = math.MaxFloat64
			}
		} else {
			w = math.MaxFloat64
		}
		z = u1 * u1 * u2
		r = gamma*v - 1.3862944
		s = a + r - w
		if s+2.609438 >= 5.0*z {
			return true
		}
		t = math.Log(z)
		return s > t
	}
	if !complete() {
		for r+alpha*math.Log(alpha/(b+w)) < t {
			if complete() {
				break
			}
		}
	}
	if aa != a {
		return b / (b + w)
	}
	return w / (b + w)
}
genBetaChengBB generates a random variate from a Beta distribution
with shape parameters aa and bb, according to Cheng's BB algorithm,
described in "Generating beta variates with non-integral shape
parameters" (1978).
package godist
import (
	"fmt"
	"math"
	"math/rand"
)
type Beta struct {
	Alpha float64
	Beta  float64
}
func (beta Beta) Mean() (float64, error) {
	if ok, err := beta.valid(); !ok {
		return 0, err
	}
	return beta.Alpha / (beta.Alpha + beta.Beta), nil
}
func (beta Beta) Median() (float64, error) {
	aa, bb := beta.Alpha, beta.Beta
	if ok, err := beta.valid(); !ok {
		return 0, err
	}
	if aa == 1 && bb > 0 {
		return 1.0 - math.Pow(2.0, -1/bb), nil
	} else if aa > 0 && bb == 1 {
		return math.Pow(2.0, -1/aa), nil
	} else if aa == 3 && bb == 2 {
		return 0.6142724318, nil
	} else if aa == 2 && bb == 3 {
		return 0.3857275681, nil
	}
	if aa == bb {
		return 0.5, nil
	} else if aa < 1 && bb < 1 {
		msg := fmt.Sprintf("Median not supported for Beta Distribution [α = %v, β = %v]", aa, bb)
		return 0, UnsupportedError{S: msg}
	}
	return (aa - 1.0/3.0) / (aa + bb - (2.0 / 3.0)), nil
}
func (beta Beta) Mode() (float64, error) {
	aa, bb := beta.Alpha, beta.Beta
	if ok, err := beta.valid(); !ok {
		return 0, err
	}
	if aa <= 1 || bb <= 1 {
		msg := fmt.Sprintf("Mode not supported for Beta Distribution [α = %v, β = %v]", aa, bb)
		return 0, InvalidDistributionError{S: msg}
	}
	return (aa - 1) / (aa + bb - 2), nil
}
func (beta Beta) Variance() (float64, error) {
	aa, bb := beta.Alpha, beta.Beta
	if ok, err := beta.valid(); !ok {
		return 0, err
	}
	return (aa * bb) / (math.Pow(aa+bb, 2) * (aa + bb + 1)), nil
}
func (beta Beta) Float64() (float64, error) {
	aa, bb := beta.Alpha, beta.Beta
	if ok, err := beta.valid(); !ok {
		return 0, err
	}
	a, b := math.Min(aa, bb), math.Max(aa, bb)
	if b < 0.5 {
		return genBetaJohnk(aa, bb), nil
	} else if a <= 1.0 {
		return genBetaChengBC(aa, bb, a, b), nil
	}
	return genBetaChengBB(aa, bb, a, b), nil
}
func genBetaJohnk(aa, bb float64) float64 {
	u, y := rand.Float64(), rand.Float64()
	return math.Pow(u, 1/aa) / (math.Pow(u, 1/aa) + math.Pow(y, 1/bb))
}