# Independently authored: direct integer binomial coefficients, no supplied imports;
# Clopper-Pearson membership by tail tests rather than endpoint bisection.
import numpy as np,math,json,sys
from pathlib import Path
p=np.arange(1,1000,dtype=float)/1000;z=1.959963984540054;s=z*z
out={m:[] for m in ['wald','wilson','agresti_coull','clopper_pearson']}
for n in range(5,101):
 k=np.arange(n+1,dtype=float);h=k/n
 mass=np.array([float(math.comb(n,int(v))) for v in k])[None,:]*p[:,None]**k[None,:]*(1-p[:,None])**(n-k)[None,:]
 center=(h+s/(2*n))/(1+s/n);half=z*np.sqrt(h*(1-h)/n+s/(4*n*n))/(1+s/n)
 a=(k+s/2)/(n+s);ah=z*np.sqrt(a*(1-a)/(n+s))
 wh=z*np.sqrt(h*(1-h)/n)
 for m,lo,hi in [('wald',h-wh,h+wh),('wilson',center-half,center+half),('agresti_coull',a-ah,a+ah)]:
  mask=(lo[None,:]<=p[:,None])&(p[:,None]<=hi[None,:]);out[m].append(np.sum(mass*mask,axis=1))
 lower=np.cumsum(mass,axis=1);upper=np.cumsum(mass[:,::-1],axis=1)[:,::-1]
 out['clopper_pearson'].append(np.sum(mass*((lower>=.025)&(upper>=.025)),axis=1))
r={}
for m,rows in out.items():
 c=np.stack(rows);r[m]={ 'below_093':round(float(np.mean(c<.93)),6),'below_095':round(float(np.mean(c<.95)),6),'mean':round(float(np.mean(c)),6),'minimum':round(float(np.min(c)),6)}
 if m=='wald':
  r[m]['fixed_p']={}
  for pv in [.01,.05,.2,.5]:
   v=c[:,round(pv*1000)-1];fall=v[:-1]-v[1:];i=int(np.argmax(fall));r[m]['fixed_p'][str(pv)]={'at_n100':round(float(v[-1]),6),'max_drop_n':i+5 if fall[i]>0 else None,'from':round(float(v[i]),6),'to':round(float(v[i+1]),6)}
print(json.dumps({'python':sys.version,'numpy':np.__version__,'independent_approximate_grid':r},indent=2))
