"""Controlled simulations, not observed Parse alert accuracy.
Re-run with Python + numpy. Fixed seed, 2,000 repetitions per scenario.
The matched-question formula follows shared/need-statistics.ts mantelHaenszelZ;
thresholds compare methods only, not all production gates or collection behavior.
"""
import csv,json,math
from pathlib import Path
import numpy as np
rng=np.random.default_rng(610062026);base=Path('public/research/ai-visibility-change-validation');N=2000;trials=[];checks=[]
def statistics(a,na,b,nb):
 total=na+nb;success=a+b
 deviation=np.sum(a-na*success/total)
 variance=np.sum(na*nb*success*(total-success)/(total*total*(total-1)))
 matched=float(deviation/np.sqrt(variance)) if variance>0 else 0.0
 diff=float(a.sum()/na.sum()-b.sum()/nb.sum());pool=float(success.sum()/total.sum())
 se=math.sqrt(pool*(1-pool)*(1/na.sum()+1/nb.sum()));pooled=diff/se if se else 0
 return diff,float(pooled),matched
for scenario in ['unchanged_independent','ten_point_increase','question_mix_only','shared_week_shocks']:
 for rep in range(N):
  p=np.linspace(.1,.7,30);na=np.full(30,14);nb=na.copy()
  if scenario=='question_mix_only':na=np.where(p>.4,14,2);nb=np.where(p>.4,2,14)
  if scenario=='shared_week_shocks':
   before_shock=rng.choice([-.1,.1],size=2);after_shock=rng.choice([-.1,.1],size=2)
   a=sum(rng.binomial(7,p+s) for s in after_shock);b=sum(rng.binomial(7,p+s) for s in before_shock)
  else:
   b=rng.binomial(nb,p);a=rng.binomial(na,p+(.1 if scenario=='ten_point_increase' else 0))
  diff,pooled,matched=statistics(a,na,b,nb)
  trials.append(dict(scenario=scenario,replicate=rep,long_run_change_points=10 if scenario=='ten_point_increase' else 0,observed_change_points=100*diff,pooled_95=int(abs(pooled)>=1.96),matched_95=int(abs(matched)>=1.96),matched_strict=int(abs(matched)>=3.29 and np.sign(diff)==np.sign(matched))))
  if rep<3:checks.append(dict(scenario=scenario,strata=[dict(current=dict(successes=int(x),trials=int(n)),previous=dict(successes=int(y),trials=int(m))) for x,n,y,m in zip(a,na,b,nb)],expected_z=matched))
metrics=[]
for scenario in sorted({r['scenario'] for r in trials}):
 rs=[r for r in trials if r['scenario']==scenario]
 for method in ['pooled_95','matched_95','matched_strict']:
  hits=sum(r[method] for r in rs);rate=hits/len(rs)
  z=1.96;n=len(rs);denom=1+z*z/n;center=(rate+z*z/(2*n))/denom;half=z*math.sqrt(rate*(1-rate)/n+z*z/(4*n*n))/denom
  metrics.append(dict(scenario=scenario,method=method,alerts=hits,trials=n,rate=rate,monte_carlo_wilson_low=max(0,center-half),monte_carlo_wilson_high=min(1,center+half)))
for name,rs in [('trials.csv',trials),('dataset.csv',metrics)]:
 with (base/name).open('w') as f:w=csv.DictWriter(f,fieldnames=list(rs[0]),lineterminator="\n");w.writeheader();w.writerows(rs)
(base/'formula-checks.json').write_text(json.dumps(checks,indent=2))
print(json.dumps({'seed':610062026,'repetitions_per_scenario':N,'questions':30,'days_per_period':14,'metrics':metrics},indent=2))