サブ　betai.js　　　Last modified: Mar 25, 2004
function betai(x, a, b)
{
  var bt

  if (x < 0.0 || x > 1.0) {
    printf("Bad x in routine betai\n")
    return 99999
  }
  if (x == 0.0 || x == 1.0) {
    bt = 0.0
  }
  else {
    bt = Math.exp(gammln(a+b)-gammln(a)-gammln(b)+a*Math.log(x)+b*Math.log(1.0-x))
  }
  if (x < (a+1.0)/(a+b+2.0)) {
    return bt*betacf(x, a, b)/a
  }
  else {
    return 1.0-bt*betacf(1.0-x, b, a)/b
  }
}


function betacf(x, a, b)
{
  var MAXIT = 100
  var EPS = 3e-7
  var FPMIN = 1e-30
  var m, m2
  var aa, c, d, del, h, qab, qam, qap

  qab = a+b
  qap = a+1
  qam = a-1
  c = 1
  d = 1-qab*x/qap
  if (Math.abs(d) < FPMIN) d = FPMIN
  d = 1/d
  h = d
  for (m = 1; m <= MAXIT; m++) {
    m2 = 2*m
    aa = m*(b-m)*x/((qam+m2)*(a+m2))
    d = 1+aa*d
    if (Math.abs(d) < FPMIN) d = FPMIN
    c = 1+aa/c
    if (Math.abs(c) < FPMIN) c = FPMIN
    d = 1/d
    h *= d*c
    aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
    d = 1+aa*d
    if (Math.abs(d) < FPMIN) d = FPMIN
    c = 1+aa/c
    if (Math.abs(c) < FPMIN) c = FPMIN
    d = 1/d
    del = d*c
    h *= del
    if (Math.abs(del-1) < EPS) break
  }
  if (m > MAXIT) {
    printf("a or b too big, or MAXIT too small in betacf\n")
    return 99999
  }
  else {
    return h
  }
}

function gammln(xx)
{
  var x, y, tmp, ser, j
  var cof = new Array(76.18009172947146, -86.50532032941677,
    24.01409824083091, -1.231739572450155,
    0.1208650973866179e-2, -0.5395239384953e-5)

  y = x = xx
  tmp = x+5.5
  tmp -= (x+0.5)*Math.log(tmp)
  ser = 1.000000000190015
  for (j = 0; j <= 5; j++) {
    ser += cof[j]/++y
  }
  return -tmp+Math.log(2.5066282746310005*ser/x)
}

function gammp(x, a)
{

  if (x < 0.0 || a <= 0.0) {
    printf("Invalid arguments in routine gammp")
    return 99999
  }
  else if (x < (a+1.0)) {
    return gser(a, x)
  }
  else {
    return 1-gcf(a, x)
  }
}


function gcf(a, x)
{
  var i, an, b, c, d, del, h, gln
  var MAXIT, EPS, FPMIN

  MAXIT = 100
  EPS = 3.0e-7
  FPMIN = 1.0e-30

  gln = gammln(a)
  b = x+1.0-a
  c = 1.0/FPMIN
  d = 1.0/b
  h = d
  for (i = 1; i <= MAXIT; i++) {
    an = -i*(i-a)
    b += 2.0
    d = an*d+b
    if (Math.abs(d) < FPMIN) d = FPMIN
    c = b+an/c
    if (Math.abs(c) < FPMIN) c = FPMIN
    d = 1.0/d
    del = d*c
    h *= del
    if (Math.abs(del-1.0) < EPS) break
  }
  if (i > MAXIT) printf("a too large, MAXIT too small in gcf")
  return Math.exp(-x+a*Math.log(x)-gln)*h
}

function gser(a, x)
{
  var n
  var sum, del, ap, gln
  var MAXIT, EPS

  MAXIT = 100
  EPS = 3.0e-7

  gln = gammln(a)
  if (x <= 0.0) {
    if (x < 0.0) printf("x less than 0 in routine gser")
    return 99999
  }
  else {
    ap = a
    del = sum = 1.0/a
    for (n = 1; n <= MAXIT; n++) {
      ++ap
      del *= x/ap
      sum += del
      if (Math.abs(del) < Math.abs(sum)*EPS) {
        return sum*Math.exp(-x+a*Math.log(x)-gln)
      }
    }
    printf("a too large, MAXIT too small in routine gser")
    return 99999
  }
}