正則化不完全ベータ関数 regularized incomplete beta function

beta_cdf.py

import math
import numpy as np
from scipy.stats import beta
import sys
args = sys.argv

# 正則化不完全ベータ関数 regularized incomplete beta function
# 0 <= z <= 1, int(a) > 0, int(b) > 0
# beta.cdf(0,a,b) == 0
# beta.cdf(1,a,b) == 1
# beta.cdf(z,a,b) == 1 - beta.cdf(1-z,b,a)

def beta_cdf1(z, a, b):
    zd = 1 - z
    t = 1
    for i in range(1, b):
        t = 1 + t * zd * (a + b - 1 - i) / (b - i)
    return math.pow(z, a) * t

def beta_cdf2(z, a, b):
    zd = (1 - z) / z
    t = 1
    for i in range(1, b):
        t = 1 + t * zd * (a + i) / (b - i)
    return math.pow(z, a + b - 1) * t

def beta_cdf3(z, a, b):
    t = 1
    for i in range(1, a):
        t = 1 + t * z * (a + b - 1 - i) / (a - i)
    return 1 - math.pow(1 - z, b) * t

def beta_cdf4(z, a, b):
    zd = z / (1 - z)
    t = 1
    for i in range(1, a):
        t = 1 + t * zd * (b + i) / (a - i)
    return 1 - math.pow(1 - z, a + b - 1) * t

z, a, b = (float(args[1]), int(args[2]), int(args[3]))
print(beta_cdf1(z, a, b))
print(beta_cdf2(z, a, b))
print(beta_cdf3(z, a, b))
print(beta_cdf4(z, a, b))
print(beta.cdf(z, a, b))