JS bayes ab calc

code.js

var g = 7;
var p = [
    0.99999999999980993,
    676.5203681218851,
    -1259.1392167224028,
    771.32342877765313,
    -176.61502916214059,
    12.507343278686905,
    -0.13857109526572012,
    9.9843695780195716e-6,
    1.5056327351493116e-7
];

var g_ln = 607/128;
var p_ln = [
    0.99999999999999709182,
    57.156235665862923517,
    -59.597960355475491248,
    14.136097974741747174,
    -0.49191381609762019978,
    0.33994649984811888699e-4,
    0.46523628927048575665e-4,
    -0.98374475304879564677e-4,
    0.15808870322491248884e-3,
    -0.21026444172410488319e-3,
    0.21743961811521264320e-3,
    -0.16431810653676389022e-3,
    0.84418223983852743293e-4,
    -0.26190838401581408670e-4,
    0.36899182659531622704e-5
];

// Spouge approximation (suitable for large arguments)
function lngamma(z) {
    if(z < 0) return Number('0/0');
    var x = p_ln[0];
    for(var i = p_ln.length - 1; i > 0; --i) x += p_ln[i] / (z + i);
    var t = z + g_ln + 0.5;
    return .5*Math.log(2*Math.PI)+(z+.5)*Math.log(t)-t+Math.log(x)-Math.log(z);
}

function lnbeta(a, b) {
  return lngamma(a) + lngamma(b) - lngamma(a+b)
  //return Math.log(gamma(a)) + Math.log(gamma(b)) + Math.log(gamma(a+b))
}

function prob_b_gt_a(alpha_a, beta_a, alpha_b, beta_b) {
  total = 0.0
  for (var i = 0; i < (alpha_b); i++) {
    total += Math.exp(lnbeta(alpha_a+i, beta_b+beta_a) - Math.log(beta_b+i) - lnbeta(1+i, beta_b) - lnbeta(alpha_a, beta_a))
  }
  return total
}