4 ms·
Quite fast and accurate: from scipy.stats import beta import numpy as np def betamax(params0, params1, integration_points=5000): u0 = para
by eutectic 6y ago
Quite fast and accurate:
from scipy.stats import beta
import numpy as np
def betamax(params0, params1, integration_points=5000):
u0 = params0[0] / sum(params0)
u1 = params1[0] / sum(params1)
flip = u0 < u1
if flip:
params0, params1 = params1, params0
q = beta(*params0).ppf(np.linspace(0, 1, integration_points + 2)[1:-1])
c = beta(*params1).cdf(q)
p = c.mean()
if flip:
return 1 - p
return p
- acidbaseextract 6y agoI appreciate the clarity in code — I'm asking for `betamaxinv` I think? >>> betamax((3, 14), (4, 12)) 0.29448379324802676 >>> betamaxinv(0.29, allowed_error=0.01) [((3, 14), (4, 12)), ...] I assume the results of `betamaxinv` are not unique, and so perhaps there need to be constraints, or return a family of results. I don't know. I'm trying to wrap my head around how OP got A ~ Beta(2,13) and B ~ Beta(3,11) without just brute forcing numbers into `betamax`.