"""Small binomial Wilson coverage versus equal-tailed Beta posterior mass; stdlib only."""
import math,json
from statistics import NormalDist
def beta_integer_cdf(x,a,b):
    n=a+b-1;return sum(math.comb(n,j)*x**j*(1-x)**(n-j) for j in range(a,n+1))
def quantile(q,a,b):
    lo,hi=0.0,1.0
    for _ in range(60):
        mid=(lo+hi)/2
        if beta_integer_cdf(mid,a,b)<q:lo=mid
        else:hi=mid
    return (lo+hi)/2
def credible(k,n,prior_a,prior_b):
    a,b=prior_a+k,prior_b+n-k;lo=quantile(.1,a,b);hi=quantile(.9,a,b);mass=beta_integer_cdf(hi,a,b)-beta_integer_cdf(lo,a,b)
    assert abs(mass-.8)<1e-12
    return {'prior_beta':[prior_a,prior_b],'posterior_beta':[a,b],'lower':lo,'upper':hi,'posterior_mass':mass,'equal_tails':[beta_integer_cdf(lo,a,b),1-beta_integer_cdf(hi,a,b)]}
def wilson(k,n):
    z=NormalDist().inv_cdf(.9);p=k/n;den=1+z*z/n;center=(p+z*z/(2*n))/den;half=z*math.sqrt(p*(1-p)/n+z*z/(4*n*n))/den
    return max(0.0,center-half),min(1.0,center+half)
def demo():
    n,k=4,1;intervals=[{'successes':j,'interval':wilson(j,n)} for j in range(n+1)];coverage=[]
    for p in [.1,.25,.5,.75,.9]:
        terms=[{'successes':j,'probability':math.comb(n,j)*p**j*(1-p)**(n-j),'contains_fixed_p':lo-1e-14<=p<=hi+1e-14} for j,(lo,hi) in enumerate(wilson(j,n) for j in range(n+1))]
        assert abs(sum(r['probability'] for r in terms)-1)<1e-12
        coverage.append({'fixed_true_p':p,'terms':terms,'actual_coverage':sum(r['probability'] for r in terms if r['contains_fixed_p'])})
    flat=credible(k,n,1,1);other=credible(k,n,2,2)
    assert flat['posterior_beta']==[2,4] and other['posterior_beta']==[3,5]
    assert beta_integer_cdf(0,2,4)==0 and beta_integer_cdf(1,2,4)==1
    assert flat['lower']!=other['lower'] and flat['upper']!=other['upper']
    return {'experiment':'four iid fictional Bernoulli trials, one observed success','n':n,'k':k,'nominal_level':.8,'wilson_observed_interval':wilson(k,n),'bayesian_flat_prior':flat,'bayesian_prior_sensitivity':other,'all_wilson_intervals':intervals,'exact_binomial_sum_coverage_no_simulation':coverage,'quantile_method':'60 bisections of exact finite polynomial integer-Beta CDF; mass tolerance 1e-12','astropy_or_scipy_executed':False}
if __name__=='__main__':print(json.dumps(demo(),ensure_ascii=False,indent=2))