Review index · Download raw source
modules/statistics_review.py
1"""Independent, lightweight count-statistics review. No flight/source-depth validation."""2from pathlib import Path3import sys, math, json, csv, ast, hashlib, datetime, argparse, shutil, subprocess4HERE=Path(__file__).resolve().parent5PACKAGE=HERE.parent6if (PACKAGE/"runtime/deps").is_dir():7 sys.path.insert(0,str(PACKAGE/"runtime/deps"))8sys.path.insert(0,str(HERE))9from common import ROOT,H,C,PHOTON_ENERGY_337,PHOTON_ENERGY_77710import numpy as np11import scipy12from scipy.stats import poisson,binom,chi2,beta13from scipy.special import xlogy,gammaln14from scipy.optimize import brentq,minimize15T_PULSE=.00316RISES_US=(30,80,200)17ALPHA=2.866515718791933e-718MARGIN=math.log(19.)1920def _integer(n):21 if isinstance(n,(bool,np.bool_)) or int(n)!=n or n<0:raise ValueError("Nonnegative integer count required")22 return int(n)23def _finite(x,nonnegative=False):24 x=float(x)25 if not math.isfinite(x) or (nonnegative and x<0):raise ValueError("Finite nonnegative value required")26 return x27def _alpha(a):28 a=_finite(a)29 if not 0<a<1:raise ValueError("0<alpha<1 required")30 return a3132def log_poisson_mass(k,mu):33 """Stable log-PMF; Stirling for k>=50. Not formal interval arithmetic."""34 k=_integer(k);mu=_finite(mu,True)35 if mu==0:return 0. if k==0 else -math.inf36 if k==0:return -mu37 if k<50:return -mu+k*math.log(mu)-math.lgamma(k+1)38 delta=(k-mu)/mu39 return (k-mu)-k*math.log1p(delta)-.5*math.log(2*math.pi*k)-1/(12*k)+1/(360*k**3)-1/(1260*k**5)4041def direct_poisson_tail(k,mu,tol=2e-14):42 """P(N>=k), positive recurrence plus decreasing geometric remainder.43 For k<=mu sum the lower tail and subtract; k1 uses stable expm1.44 Includes no Gaussian approximation. PMF roundoff is separate from tail truncation.45 """46 k=_integer(k);mu=_finite(mu,True)47 if k==0:return 1.48 if mu==0:return 0.49 if k==1:return -math.expm1(-mu)50 if k-1<mu:51 j=k-1;term=math.exp(log_poisson_mass(j,mu));terms=[term];s=term52 while j>0:53 ratio=j/mu54 if ratio<1 and term*ratio/(1-ratio)<=tol*s:break55 term*=ratio;j-=1;terms.append(term);s+=term56 return max(0.,min(1.,1-math.fsum(terms)))57 j=k;term=math.exp(log_poisson_mass(j,mu));terms=[term];s=term58 while True:59 ratio=mu/(j+1)60 if term==0 or term*ratio/(1-ratio)<=tol*s:break61 j+=1;term*=mu/j;terms.append(term);s+=term62 if len(terms)>200000:raise RuntimeError("Use a qualified large-parameter method")63 return max(0.,min(1.,math.fsum(terms)))6465def critical_known(b,alpha=ALPHA,method="scipy"):66 """Smallest integer k>=1 satisfying specified single-test tail."""67 b=_finite(b,True);alpha=_alpha(alpha)68 tail=(lambda k:float(poisson.sf(k-1,b))) if method=="scipy" else (lambda k:direct_poisson_tail(k,b))69 lo=0;hi=max(1,math.ceil(b+12*math.sqrt(b+1)+60))70 while tail(hi)>alpha:hi*=271 while hi-lo>1:72 m=(hi+lo)//273 if tail(m)<=alpha:hi=m74 else:lo=m75 return hi7677def required_signal_known(b,alpha=ALPHA,power=.95):78 """Expected detected source count s; known b in same ON window."""79 b=_finite(b,True)80 if not 0<power<1:raise ValueError("0<power<1 required")81 k=critical_known(b,alpha)82 hi=max(b+10,k+10*math.sqrt(k+1))83 while poisson.sf(k-1,hi)<power:hi*=284 total=b if poisson.sf(k-1,b)>=power else brentq(lambda m:float(poisson.sf(k-1,m))-power,b,hi,xtol=1e-11,rtol=1e-14)85 return {"critical":k,"source_mean_PE":total-b,"background_mean_PE":b,86 "power":float(poisson.sf(k-1,total)),"tail":float(poisson.sf(k-1,b))}8788def onoff_pvalue(n,m,r):89 """N~Pois(s+b), M~Pois(r*b); r=T_off/T_on."""90 n=_integer(n);m=_integer(m);r=_finite(r)91 if r<=0:raise ValueError("r=T_off/T_on must be positive")92 return float(binom.sf(n-1,n+m,1/(1+r)))9394def direct_onoff_tail(n,m,r):95 n=_integer(n);m=_integer(m);r=_finite(r)96 if r<=0:raise ValueError("Positive exposure ratio required")97 t=n+m;p=1/(1+r)98 if n==0:return 1.99 terms=[math.exp(math.lgamma(t+1)-math.lgamma(j+1)-math.lgamma(t-j+1)100 +j*math.log(p)+(t-j)*math.log1p(-p)) for j in range(n,t+1)]101 return math.fsum(terms)102103def source_upper_known(n,b,alpha=.05):104 """At-least nominal unconditional coverage for expected source PE."""105 return max(0.,float(chi2.isf(_alpha(alpha),2*(_integer(n)+1))/2)-_finite(b,True))106107def source_upper_off(n,m,r,alpha=.05,split=.5):108 """Conservative union-bound UL, fixed alpha split; r=T_off/T_on."""109 n=_integer(n);m=_integer(m);r=_finite(r);alpha=_alpha(alpha);split=_alpha(split)110 if r<=0:raise ValueError("Positive exposure ratio required")111 U=float(chi2.isf(alpha*split,2*(n+1))/2)112 L=0. if m==0 else float(chi2.ppf(alpha*(1-split),2*m)/2)113 return max(0.,U-L/r)114115def mc_intervals(k,n,alpha=.05):116 """Individual Wilson and exact Clopper-Pearson binomial intervals."""117 k=_integer(k);n=_integer(n);alpha=_alpha(alpha)118 if k>n:raise ValueError("k<=n required")119 if n==0:return {"fraction":None,"Wilson":[None,None],"Clopper_Pearson":[None,None]}120 from scipy.stats import norm121 z=float(norm.isf(alpha/2));p=k/n;den=1+z*z/n122 mid=(p+z*z/(2*n))/den;half=z*math.sqrt(p*(1-p)/n+z*z/(4*n*n))/den123 cp=[0. if k==0 else float(beta.ppf(alpha/2,k,n-k+1)),124 1. if k==n else float(beta.isf(alpha/2,k+1,n-k))]125 return {"fraction":p,"Wilson":[max(0.,mid-half),min(1.,mid+half)],"Clopper_Pearson":cp}126127def primitive(t,rise_us):128 z=np.clip(np.asarray(t,float),0,T_PULSE);d=rise_us*1e-6/.8129 return np.where(z<d,z*z/(2*d),z-d/2)130131def q_integrals(left,right,onset,rise_us):132 return primitive(np.asarray(right)-onset,rise_us)-primitive(np.asarray(left)-onset,rise_us)133134def poisson_loglik(n,mu):135 n=np.asarray(n,float);mu=np.asarray(mu,float)136 if np.any(mu<0):return -math.inf137 return float(np.sum(xlogy(n,mu)-mu-gammaln(n+1)))138139def nuisance_profile(n,d,q):140 """Independent full-bin concave mixture solution, not imported R4 code."""141 n=np.asarray(n,float);d=np.asarray(d,float);q=np.asarray(q,float)142 if np.any(n<0) or np.any(d<=0) or np.any(q<0):raise ValueError("Invalid counts/windows")143 N=float(n.sum());D=float(d.sum());Q=float(q.sum())144 if Q<=0:raise ValueError("Nonzero source support required")145 if N==0:return {"theta":0.,"A_cps":0.,"b_cps":0.,"logL":0.,"mu_hat":np.zeros(len(n))}146 h=q/Q;w=d/D;a=h/w-1147 def derivative(x):148 use=n>0149 return float(np.sum(n[use]*a[use]/(1+x*a[use])))150 if derivative(0)<=0:theta=0.151 elif np.all(n[q==0]==0) and derivative(1-1e-14)>=0:theta=1.152 else:theta=brentq(derivative,0,1-1e-14,xtol=1e-14,rtol=1e-14)153 A=theta*N/Q;b=(1-theta)*N/D;mu=A*q+b*d154 return {"theta":theta,"A_cps":A,"b_cps":b,"logL":poisson_loglik(n,mu),"mu_hat":mu}155156def worked_example(seed=202610221110):157 """Actual generated trial with explicit n,d,q and an independent 2D fit."""158 left=np.r_[np.arange(-100,400,20)*1e-6,.001,-.0105]159 right=np.r_[np.arange(-80,420,20)*1e-6,.0025,-.0005]160 d=right-left;truth={"rise_us":80,"onset_s":35e-6,"A_cps":1000/.003,"b_cps":1e4}161 q=q_integrals(left,right,truth["onset_s"],truth["rise_us"])162 q[-1]=0.163 rng=np.random.default_rng(seed);mu=truth["A_cps"]*q+truth["b_cps"]*d164 n=rng.poisson(mu);fit=nuisance_profile(n,d,q)165 # Direct original lambda=A*q+b*d optimized in two independent rate coordinates.166 scales=np.array([1e5,1e4])167 def obj(x):return -poisson_loglik(n,(x[0]*scales[0])*q+(x[1]*scales[1])*d)168 def jac(x):169 m=(x[0]*scales[0])*q+(x[1]*scales[1])*d170 ratio=np.divide(n,m,out=np.zeros_like(m),where=m>0)-1171 return -np.array([np.sum(ratio*q)*scales[0],np.sum(ratio*d)*scales[1]])172 opt=minimize(obj,[2.,1.],jac=jac,method="L-BFGS-B",bounds=[(1e-12,None)]*2,173 options={"ftol":1e-15,"gtol":1e-9,"maxiter":300})174 direct_rates=opt.x*scales175 assert abs(fit["logL"]+opt.fun)<2e-7176 assert np.max(abs(direct_rates/np.array([fit["A_cps"],fit["b_cps"]])-1))<2e-5177 rows=[]178 for i in range(len(n)):179 rows.append({"bin":i,"label":"prewindow" if i==len(n)-1 else "late_plateau" if i==len(n)-2 else "early",180 "left_s":float(left[i]),"right_s":float(right[i]),"d_s":float(d[i]),"q_s":float(q[i]),181 "n_PE":int(n[i]),"truth_mean_PE":float(mu[i]),"fitted_mean_PE":float(fit["mu_hat"][i])})182 candidates=[]183 for rise in RISES_US:184 for onset in np.arange(-50,100.1,5)*1e-6:185 qq=q_integrals(left,right,onset,rise);qq[-1]=0186 f=nuisance_profile(n,d,qq)187 candidates.append({"rise_us":rise,"onset_us":float(onset*1e6),"A_cps":f["A_cps"],"b_cps":f["b_cps"],"logL":f["logL"]})188 candidates.sort(key=lambda x:x["logL"],reverse=True)189 result={"status":"EXECUTED_ONE_SYNTHETIC_TRIAL_NOT_EMPIRICAL_SOURCE_DATA","seed":seed,"truth":truth,190 "observed_total_PE":int(n.sum()),"window_total_s":float(d.sum()),"source_integral_q_s":float(q.sum()),191 "profile":{"theta":fit["theta"],"A_cps":fit["A_cps"],"b_cps":fit["b_cps"],"logL":fit["logL"]},192 "direct_SciPy_two_parameter_fit":{"success":bool(opt.success),"message":str(opt.message),"A_cps":float(direct_rates[0]),"b_cps":float(direct_rates[1]),"logL":float(-opt.fun)},193 "absolute_loglikelihood_difference":abs(fit["logL"]+opt.fun),194 "fitted_total_minus_observed_PE":float(fit["mu_hat"].sum()-n.sum()),195 "bins":rows,"top_grid_candidates":candidates[:9],196 "qualification":"The true80us candidate is used for the explicit2D/profile identity check. The separate3rise/onset search illustrates unknown morphology and onset; fitted labels are not source depth."}197 return result198199def photon_planes(source_PE,efficiency=.2,diameter_mm=5):200 """Aperture fluence link; expected primary PE before multiplication."""201 eta=_finite(efficiency);s=_finite(source_PE,True);d=_finite(diameter_mm)202 if not 0<eta<=1 or d<=0:raise ValueError("Positive aperture/efficiency required")203 A=math.pi*(d*1e-3)**2/4204 return {"expected_source_primary_PE":s,"assumed_efficiency":eta,"area_m2":A,205 "expected_aperture_photons":s/eta,"aperture_photon_fluence_m2":s/(eta*A),206 "aperture_energy_fluence_J_m2_at337":s/(eta*A)*PHOTON_ENERGY_337,207 "E337_J":PHOTON_ENERGY_337,"E777_J":PHOTON_ENERGY_777,208 "qualification":"Aperture photons after propagation; no emitted-source yield or second cloud attenuation; no PMT/APD gain factor."}209210def count_checks():211 rows=[]212 for mu in [0.,1e-12,.01,.1,1.,10.,100.,1000.,1e5,1e6]:213 ns=sorted({0,1,max(0,int(mu)),int(mu+5*math.sqrt(mu+1))})214 for n in ns:215 a=direct_poisson_tail(n,mu);b=float(poisson.sf(n-1,mu))216 assert abs(a-b)<1e-8217 rows.append({"kind":"Poisson_tail","mean":mu,"n_inclusive":n,"direct_sum":a,"SciPy":b,218 "abs_difference":abs(a-b),"relative_difference":abs(a/b-1) if b else None})219 thresholds=[]220 for b in [0.,.01,1.,10.,100.,1000.,1e5,1e6]:221 for a in [ALPHA,.1/(50000*86400*3)]:222 t=required_signal_known(b,a);k2=critical_known(b,a,method="sum")223 assert k2==t["critical"]224 direct_power=direct_poisson_tail(k2,b+t["source_mean_PE"])225 assert abs(direct_power-.95)<1e-8226 thresholds.append({"background_PE":b,"alpha_per_test":a,**t,227 "critical_direct_sum":k2,"power_direct_sum":direct_power})228 onoff=[]229 for n,m,r in [(0,0,10),(1,0,10),(10,0,10),(5,2,3),(30,40,10),(100,200,10)]:230 a=onoff_pvalue(n,m,r);b=direct_onoff_tail(n,m,r);assert abs(a-b)<1e-12231 onoff.append({"n_on":n,"m_off":m,"r_off_over_on":r,"SciPy":a,"finite_sum":b})232 assert abs(source_upper_known(0,0)+math.log(.05))<1e-13233 assert abs(source_upper_off(0,0,10)+math.log(.025))<1e-13234 upper=[]235 for n,b in [(0,0),(0,2),(1,0),(3,2),(10,5),(100,80)]:236 u=source_upper_known(n,b)237 total=float(chi2.isf(.05,2*(n+1))/2)238 cdf_direct=1-direct_poisson_tail(n+1,total)239 assert abs(cdf_direct-.05)<1e-11240 upper.append({"n_on":n,"known_b_PE":b,"source95_upper_PE":u,"total_mean_inversion_CDF":cdf_direct})241 # Coverage is an independent finite enumeration, not a Monte Carlo claim.242 coverage=[]243 for src,b,r in [(0,0,10),(1,0,10),(3,2,10),(10,3,5)]:244 n=np.arange(100);m=np.arange(200)245 pn=poisson.pmf(n,src+b);pm=poisson.pmf(m,r*b)246 UK=np.maximum(0.,chi2.isf(.05,2*(n+1))/2-b)247 up_total=chi2.isf(.025,2*(n+1))/2248 lower=np.zeros(len(m));lower[1:]=chi2.ppf(.025,2*m[1:])/2249 UF=np.maximum(0.,up_total[:,None]-lower[None,:]/r)250 pk=float(pn[UK>=src-1e-12].sum())251 pf=float(np.sum((UF>=src-1e-12)*pn[:,None]*pm[None,:]))252 omitted=max(0.,1-float(pn.sum()*pm.sum()))253 assert pk>=.95-1e-10 and pf>=.95-1e-10254 coverage.append({"source_PE":src,"b_on_PE":b,"r":r,"known_b_enumerated_coverage":pk,255 "finite_off_enumerated_coverage":pf,"omitted_probability_bound":omitted})256 pins=[]257 for mu,k,pinned in [(1000000.,1005000,2.934034048031641e-7),258 (1000000.,1007000,1.3595737503819116e-12),259 (1000000.,1005004,2.8739183629417515e-7),260 (1000000.,1005005,2.859076013349398e-7)]:261 val=direct_poisson_tail(k,mu);rel=abs(val/pinned-1)262 assert rel<1e-8263 pins.append({"mean":mu,"n_inclusive":k,"pinned_60digit_Decimal_tail":pinned,264 "direct_recurrence":val,"relative_difference":rel,265 "SciPy_tail":float(poisson.sf(k-1,mu)),266 "provenance":"R4 qa/browser_large_mean_tail_control.json;60digit Decimal/Stirling/geometric-remainder calculation, not empirical data."})267 return {"tail_checks":rows,"pinned_high_mean_controls":pins,"threshold_power_checks":thresholds,"onoff_checks":onoff,268 "upper_inversions":upper,"coverage_finite_enumerations":coverage,269 "zero_count_controls":{"known_zero_background95":source_upper_known(0,0),"finite_baseline_zero_zero95":source_upper_off(0,0,10)},270 "qualification":"Exact probability constructions, floating-point implementations. At mean1e6 the SciPy incomplete-gamma tail differs from direct positive summation by several parts per million near5sigma; integer checks here agree. This is not formal interval arithmetic."}271272def confusion_checks():273 report=json.loads((ROOT/"data/statistics_results.json").read_text(encoding="utf-8-sig"))274 pr=report["practical_morphology"];rows=[]275 metric_count=0;max_ci_difference=0.276 for case in pr["scenarios"]:277 trial=ROOT/"data/trials"/f"case_{case['case_id']:03d}.npz"278 if not trial.is_file():raise FileNotFoundError("Missing synthetic trial summaries: "+str(trial))279 with np.load(trial,allow_pickle=False) as t:280 scores=t["scores"];truth=RISES_US.index(case["truth_receiver_rise_10_90_us"])281 best=scores.argmax(axis=1);group=(scores[:,2]>scores[:,:2].max(axis=1)).astype(int)282 detected=np.asarray(t["detection_pvalue"])<=ALPHA283 pcheck=binom.sf(t["gate_count"]-1,t["gate_count"]+t["pre_count"],.0034/(.0034+.010))284 assert np.allclose(pcheck,t["detection_pvalue"],rtol=1e-8,atol=1e-14)285 sort=np.sort(scores,axis=1);margin3=sort[:,-1]-sort[:,-2]286 margin2=np.abs(scores[:,2]-scores[:,:2].max(axis=1))287 a3=detected&(margin3>=MARGIN);a2=detected&(margin2>=MARGIN)288 pred3=np.where(~detected,4,np.where(~a3,3,best))289 pred2=np.where(~detected,3,np.where(~a2,2,group))290 c3=np.bincount(pred3,minlength=5).tolist();c2=np.bincount(pred2,minlength=4).tolist()291 assert c3==case["three_class"]["confusion_counts"] and c2==case["fast_slow"]["confusion_counts"]292 assert len(best)==case["trials"]==1000293 metrics=[(case["detection"],int(detected.sum()),len(best)),294 (case["three_class"]["accepted"],int(a3.sum()),len(best)),295 (case["three_class"]["wrong_among_accepted"],int((a3&(best!=truth)).sum()),int(a3.sum())),296 (case["fast_slow"]["accepted"],int(a2.sum()),len(best)),297 (case["fast_slow"]["wrong_among_accepted"],int((a2&(group!=(truth==2))).sum()),int(a2.sum())),298 (case["three_class"]["rejected_among_detected"],int((detected&~a3).sum()),int(detected.sum())),299 (case["fast_slow"]["rejected_among_detected"],int((detected&~a2).sum()),int(detected.sum())),300 (case["unthresholded_practical_three_class_error"],int((best!=truth).sum()),len(best))]301 rebuilt=[]302 for old,k,n in metrics:303 assert old["count"]==k and old["denominator"]==n304 ii=mc_intervals(k,n);metric_count+=1305 if n:306 for field,newkey in [("individual_95pct_Wilson_interval","Wilson"),307 ("individual_95pct_Clopper_Pearson_interval","Clopper_Pearson")]:308 if field in old:309 diff=float(np.max(abs(np.array(old[field])-ii[newkey])))310 max_ci_difference=max(max_ci_difference,diff);assert diff<2e-12311 rebuilt.append({"k":k,"n":n,**ii})312 rows.append({"case_id":case["case_id"],"N_plateau_equivalent_PE":case["N_plateau_equivalent_3ms_PE"],313 "background_cps":case["background_cps"],"sample_khz":case["sample_rate_khz"],314 "truth_rise_us":case["truth_receiver_rise_10_90_us"],"confusion3":c3,"confusion_fastslow":c2,315 "metrics_order":["detected","accepted3","wrong_given_accepted3","accepted_fastslow","wrong_given_accepted_fastslow","rejected3_given_detected","rejected_fastslow_given_detected","unthresholded3_error"],316 "rebuilt_metrics":rebuilt,"trial_sha256":hashlib.sha256(trial.read_bytes()).hexdigest()})317 assert len(rows)==81318 return {"status":"ALL_SAVED_TRIAL_SUMMARIES_RECOMPUTED","scenarios":81,"signal_trials":81000,319 "metrics_recomputed":metric_count,"max_interval_absolute_difference":max_ci_difference,"rows":rows,320 "qualification":"Rebuilds decisions and denominators from saved synthetic fitted scores, not raw experimental traces. Does not independently refit every pulse because full per-bin count arrays were not retained in these summaries. No depth-label or flight validation."}321322def original_source_map():323 requested={"oracle_reference/calculate_ideal_poisson_rise.py":["experiment","affinity_bounds"],"practical_morphology.py":["primitive","fit_one","draw_case","analyze_case"],324 "count_inference/count_inference.py":["count_threshold_known_background","exact_onoff_pvalue",325 "poisson_total_mean_upper","poisson_total_mean_lower","red_source_upper_known_background",326 "red_source_upper_finite_baseline"]}327 result=[]328 base=ROOT/"reference/statistics"329 for rel,names in requested.items():330 p=base/rel331 if not p.is_file():result.append({"file":str(p.relative_to(ROOT)),"present":False});continue332 tree=ast.parse(p.read_text(encoding="utf-8-sig"))333 line={n.name:n.lineno for n in ast.walk(tree) if isinstance(n,(ast.FunctionDef,ast.AsyncFunctionDef))}334 result.append({"file":str(p.relative_to(ROOT)),"present":True,"sha256":hashlib.sha256(p.read_bytes()).hexdigest(),335 "functions":{n:line.get(n) for n in names}})336 return result337338def javascript_crosscheck(node_path=None):339 p=ROOT/"reference/statistics/browser_math.js"340 if not p.is_file():return {"status":"NOT_RUN","reason":"Reference browser_math.js not present"}341 lines=p.read_text(encoding="utf-8-sig").splitlines()342 mapping={name:next((i+1 for i,line in enumerate(lines) if "function "+name+"(" in line),None)343 for name in ["sf","critical","threshold","upper"]}344 node=node_path or shutil.which("node")345 receipt={"file":str(p.relative_to(ROOT)),"sha256":hashlib.sha256(p.read_bytes()).hexdigest(),"function_lines":mapping}346 if not node:return {**receipt,"status":"NOT_RUN","reason":"Node optional and unavailable; Python checks do not execute JavaScript"}347 cases=[]348 for b in [0.,.01,1.,10.,1000.,1e6]:349 for alpha in [ALPHA,.1/(50000*86400*3)]:cases.append([b,alpha])350 red_cases=[[0,0],[0,2],[1,0],[3,2],[10,5],[100,80]]351 tiny_cases=[]352 for alpha in [ALPHA,.1/(50000*86400*3)]:353 for factor in [1-1e-8,1+1e-8]:354 tiny_cases.append([-math.log1p(-alpha)*factor,alpha])355 script="const M=require(process.argv[1]);const q=JSON.parse(process.argv[2]);process.stdout.write(JSON.stringify({threshold:q.threshold.map(a=>M.threshold(...a)),upper:q.upper.map(a=>M.upper(...a)),tiny:q.tiny.map(a=>({sf1:M.sf(1,a[0]),critical:M.critical(...a)}))}));"356 payload={"threshold":cases,"upper":red_cases,"tiny":tiny_cases}357 run=subprocess.run([node,"-e",script,str(p),json.dumps(payload)],text=True,capture_output=True,check=True)358 outputs=json.loads(run.stdout);checks=[]359 for (b,a),js in zip(cases,outputs["threshold"]):360 ref=required_signal_known(b,a)361 assert js["critical"]==ref["critical"] and abs(js["signal"]-ref["source_mean_PE"])<5e-5362 checks.append({"b_PE":b,"alpha":a,"JS":js,"Python":ref,"signal_abs_difference":abs(js["signal"]-ref["source_mean_PE"])})363 red=[]364 for (n,b),js in zip(red_cases,outputs["upper"]):365 ref=source_upper_known(n,b);assert abs(js-ref)<1e-8366 red.append({"n_on":n,"b_on_PE":b,"JS_upper95":js,"Python_upper95":ref,"abs_difference":abs(js-ref)})367 tiny=[]368 for (b,a),js in zip(tiny_cases,outputs["tiny"]):369 ref=-math.expm1(-b);critical=critical_known(b,a,method="sum")370 assert js["critical"]==critical and abs(js["sf1"]/ref-1)<1e-12371 tiny.append({"b_PE":b,"alpha":a,"JS":js,"expected_sf1":ref,"expected_critical":critical})372 version=subprocess.run([node,"--version"],text=True,capture_output=True,check=True).stdout.strip()373 return {**receipt,"status":"EXECUTED_OPTIONAL_NODE","node_version":version,"checks":checks,374 "red_upper_checks":red,"tiny_background_boundary_checks":tiny}375376def profile_boundary_checks():377 rows=[]378 for n,theta in [([0,0],0.),([10,0],0.),([0,10],1.),([10,20],1/3)]:379 fit=nuisance_profile(n,[1.,1.],[0.,1.])380 assert abs(fit["theta"]-theta)<1e-12381 assert abs(float(fit["mu_hat"].sum())-sum(n))<1e-12382 rows.append({"counts":n,"expected_theta":theta,"fitted_theta":fit["theta"],383 "A_cps":fit["A_cps"],"b_cps":fit["b_cps"],"logL":fit["logL"]})384 return rows385386def affinity_bound_values(mu,nu):387 mu=np.asarray(mu,float);nu=np.asarray(nu,float)388 if mu.shape!=nu.shape or np.any(mu<0) or np.any(nu<0):raise ValueError("Matched nonnegative bin means required")389 h=float(.5*np.sum((np.sqrt(mu)-np.sqrt(nu))**2))390 bc=math.exp(-h)391 root=math.sqrt(max(0.,1-bc*bc))392 return {"affinity":bc,"Bayes_error_lower":bc*bc/(2*(1+root)),"Bayes_error_upper":bc/2}393394def oracle_affinity_controls():395 rows=[]396 for mu,nu in [([.2,1.],[1.3,2.5]),([0.,.3],[.2,1.]),([1.,.1],[1.,3.])]:397 pa=[np.array([math.exp(log_poisson_mass(k,m)) for k in range(60)]) for m in mu]398 qa=[np.array([math.exp(log_poisson_mass(k,m)) for k in range(60)]) for m in nu]399 p=np.outer(pa[0],pa[1]);q=np.outer(qa[0],qa[1])400 bc=float(np.sqrt(p*q).sum());pe=float(.5*np.minimum(p,q).sum())401 formula=affinity_bound_values(mu,nu)402 assert abs(p.sum()-1)<1e-13 and abs(q.sum()-1)<1e-13403 assert abs(bc-formula["affinity"])<1e-13404 assert formula["Bayes_error_lower"]-1e-13<=pe<=formula["Bayes_error_upper"]+1e-13405 rows.append({"mu_PE_per_bin":mu,"nu_PE_per_bin":nu,"explicit_joint_PMF_affinity":bc,406 "explicit_equal_prior_Bayes_error":pe,**formula,"finite_grid_mass_P":float(p.sum()),"finite_grid_mass_Q":float(q.sum())})407 return {"status":"THREE_NEW_TWO_BIN_CONTROLS_EXECUTED","rows":rows,408 "retained320_grid":"NOT_INDEPENDENTLY_REBUILT_IN_THIS_LIGHT_REVIEW",409 "assumptions":"Known simple hypotheses, independent Poisson bins, equal priors and0-1loss. Unknown nuisance fits or rejected/accepted-only errors do not inherit these bounds."}410411def _numeric_assignments(path):412 """Read simple numeric constants with AST; never execute retained scripts."""413 env={}414 def calc(n):415 if isinstance(n,ast.Constant) and isinstance(n.value,(int,float)):return n.value416 if isinstance(n,ast.Name):return env[n.id]417 if isinstance(n,ast.BinOp):418 x=calc(n.left);y=calc(n.right)419 if isinstance(n.op,ast.Mult):return x*y420 if isinstance(n.op,ast.Div):return x/y421 if isinstance(n.op,ast.Add):return x+y422 if isinstance(n.op,ast.Sub):return x-y423 raise ValueError("Not a simple numeric expression")424 for node in ast.parse(path.read_text(encoding="utf-8-sig")).body:425 if isinstance(node,ast.Assign) and len(node.targets)==1 and isinstance(node.targets[0],ast.Name):426 try:env[node.targets[0].id]=calc(node.value)427 except (ValueError,KeyError):pass428 return env429430def downstream_threshold_bindings():431 data=json.loads((ROOT/"data/analytic_results.json").read_text(encoding="utf-8-sig"))432 ref=ROOT/"reference/analytic/count_requirements.py"433 fovref=ROOT/"reference/analytic/fov_population.py"434 constants=_numeric_assignments(ref)435 operational=constants["ALPHA_OPERATIONAL"]436 assert abs(operational-.1/(50000*86400*3))<1e-25437 source=fovref.read_text(encoding="utf-8-sig")438 assert "m.count_at_power((rate+dark)*dt,m.ALPHA_OPERATIONAL,.95)" in source439 dark=float(data["fov_population"]["report"]["inputs"]["dark_cps"])440 cache={};rows=[]441 for study in ["fov_population","auroral_counts"]:442 for index,row in enumerate(data[study]["rows"]):443 if study=="fov_population":444 b=(row["background_photoelectron_cps"]+dark)*row["window_s"]445 alpha=operational;storedS=row["required_source_count_for_95pct_detection"]446 alpha_source="reference/analytic/fov_population.py:41 -> count_requirements.py ALPHA_OPERATIONAL"447 else:448 b=row["known_background_mean_count"];alpha=row["per_test_alpha"]449 storedS=row["signal_mean_for_95pct_power"];alpha_source="row.per_test_alpha"450 key=(b,alpha)451 if key not in cache:cache[key]=required_signal_known(b,alpha)452 check=cache[key]453 delta=abs(storedS-check["source_mean_PE"])454 assert row["critical_observed_count"]==check["critical"]455 assert delta<=5e-5+1e-9*abs(storedS)456 achieved=float(poisson.sf(check["critical"]-1,b+storedS))457 assert abs(achieved-.95)<1e-8458 rows.append({"study":study,"row_index":index,"background_mean_PE":b,459 "alpha_per_test":alpha,"alpha_source":alpha_source,460 "stored_critical":row["critical_observed_count"],"recomputed_critical":check["critical"],461 "stored_S95_PE":storedS,"recomputed_S95_PE":check["source_mean_PE"],462 "S95_abs_difference_PE":delta,"power_at_stored_S95":achieved})463 assert sum(x["study"]=="fov_population" for x in rows)==216464 assert sum(x["study"]=="auroral_counts" for x in rows)==192465 return {"status":"ALL408_STATISTICAL_COLUMN_BINDINGS_CHECKED",466 "rows_checked":len(rows),"unique_background_alpha_computations":len(cache),467 "FOV_rows":216,"auroral_rows":192,468 "maximum_S95_absolute_difference_PE":max(x["S95_abs_difference_PE"] for x in rows),469 "maximum_power_absolute_difference":max(abs(x["power_at_stored_S95"]-.95) for x in rows),470 "constant_provenance":{"file":str(ref.relative_to(ROOT)),"sha256":hashlib.sha256(ref.read_bytes()).hexdigest(),471 "ALPHA_OPERATIONAL":operational,"dark_cps_from_FOV_report":dark},472 "qualification":"Checks only known-background Poisson critical/S95columns and their binding to actualB/declared false-alarm alpha. Population-slope alpha is a different variable. Does not validate FOVsource geography, brightness law, background scene or mission rate.",473 "rows":rows}474475476477def equation_function_map():478 mapping={479 "S01":["photon_planes"],"S02":["photon_planes"],480 "S03":["onoff_pvalue"],"S04":["onoff_pvalue"],"S05":["critical_known"],481 "S06":["required_signal_known"],"S07":["onoff_pvalue","direct_onoff_tail"],482 "S08":["count_checks"],"S09":["source_upper_known"],"S10":["source_upper_known"],483 "S11":["source_upper_off"],"S12":["source_upper_off"],484 "S13":["q_integrals","poisson_loglik"],"S14":["nuisance_profile"],485 "S15":["nuisance_profile"],"S16":["nuisance_profile","worked_example"],486 "S17":["confusion_checks"],"S18":["mc_intervals"],"S19":["mc_intervals"],"S20":["affinity_bound_values"],"S21":["oracle_affinity_controls"]}487 tree=ast.parse(Path(__file__).read_text(encoding="utf-8-sig"))488 lines={n.name:n.lineno for n in ast.walk(tree) if isinstance(n,ast.FunctionDef)}489 return {eq:[{"file":"modules/statistics_review.py","function":f,"line":lines[f]} for f in funs]490 for eq,funs in mapping.items()}491492def run_checks(write=False,node_path=None):493 checks={"schema":"INDEPENDENT_STATISTICS_REVIEW_V1","executed_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),494 "versions":{"python":sys.version.split()[0],"numpy":np.__version__,"scipy":scipy.__version__},495 "photon_planes":photon_planes(100),496 "counts":count_checks(),"oracle_affinity":oracle_affinity_controls(),"downstream_threshold_bindings":downstream_threshold_bindings(),"profile_boundaries":profile_boundary_checks(),"worked_example":worked_example(),"confusions":confusion_checks(),497 "source_map":original_source_map(),"equation_function_map":equation_function_map(),"javascript":javascript_crosscheck(node_path),498 "theory_scope":"Own derivations are in methods/statistics.md. Independent numerical comparisons check implementation, not physical assumptions or source classes.",499 "mission_rate":"NOT_COMPUTABLE_FROM_THESE_INPUTS","passed":True}500 if write:501 out=ROOT/"checks";out.mkdir(exist_ok=True)502 (out/"statistics_results.json").write_text(json.dumps(checks,indent=2,allow_nan=False),encoding="utf-8")503 rows=checks["worked_example"]["bins"]504 with (out/"statistics_worked_bins.csv").open("w",newline="",encoding="utf-8") as f:505 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)506 (out/"statistics_summary.json").write_text(json.dumps({507 "passed":True,"worked_logL_abs_difference":checks["worked_example"]["absolute_loglikelihood_difference"],508 "scenarios_rebuilt":checks["confusions"]["scenarios"],"metrics_rebuilt":checks["confusions"]["metrics_recomputed"],509 "max_interval_abs_difference":checks["confusions"]["max_interval_absolute_difference"],510 "zero_count_controls":checks["counts"]["zero_count_controls"],"javascript_status":checks["javascript"]["status"]},indent=2),encoding="utf-8")511 return checks512513if __name__=="__main__":514 arg=argparse.ArgumentParser();arg.add_argument("--write",action="store_true");arg.add_argument("--node")515 a=arg.parse_args();result=run_checks(a.write,a.node)516 print(json.dumps({"passed":result["passed"],"worked_logL_difference":result["worked_example"]["absolute_loglikelihood_difference"],517 "scenarios":result["confusions"]["scenarios"],"zero_counts":result["counts"]["zero_count_controls"],518 "JavaScript":result["javascript"]["status"]},indent=2))519520