fn ln_beta(a: f64, b: f64) -> f64
Log of the Beta function: ln(B(a,b)) = ln(Gamma(a)) + ln(Gamma(b)) - ln(Gamma(a+b))