"""Illustrative short-form precision model; Python standard library only. No fitting or item selection is performed. See the simulator page for limitations. Power and required N use the exact noncentral t distribution (two-sided alpha = .05, two equal groups), matching R's pwr.t.test. """ import functools import math import json import sys from statistics import NormalDist def _cdf(x): return .5*math.erfc(-x/math.sqrt(2)) def _betacf(x, a, b): tiny = 1e-300 c, d = 1.0, 1 - (a+b)*x/(a+1) d = 1/(d if abs(d) >= tiny else tiny) h = d for m in range(1, 501): m2 = 2*m for aa in (m*(b-m)*x/((a+m2-1)*(a+m2)), -(a+m)*(a+b+m)*x/((a+m2)*(a+m2+1))): d = 1 + aa*d d = 1/(d if abs(d) >= tiny else tiny) c = 1 + aa/c c = c if abs(c) >= tiny else tiny delta = d*c h *= delta if abs(delta-1) < 1e-15: break return h def _ibeta(x, a, b): """Regularized incomplete beta I_x(a, b).""" if x <= 0: return 0.0 if x >= 1: return 1.0 bt = math.exp(math.lgamma(a+b)-math.lgamma(a)-math.lgamma(b)+a*math.log(x)+b*math.log(1-x)) if x < (a+1)/(a+b+2): return bt*_betacf(x, a, b)/a return 1 - bt*_betacf(1-x, b, a)/b def _pt(t, df): tail = .5*_ibeta(df/(df+t*t), df/2, .5) return 1 - tail if t >= 0 else tail def _dt(x, df): return math.exp(math.lgamma((df+1)/2)-math.lgamma(df/2)-.5*math.log(df*math.pi) -(df+1)/2*math.log1p(x*x/df)) @functools.lru_cache(maxsize=None) def _qt(p, df): """Upper quantile (p > .5): Newton's method from a Cornish-Fisher start, bisection fallback.""" z = NormalDist().inv_cdf(p) x, ok = z + (z**3+z)/(4*df), False for _ in range(50): if not (x > 0 and math.isfinite(x)): break dx = (_pt(x, df)-p)/_dt(x, df) x -= dx if abs(dx) < 1e-12*max(1, abs(x)): ok = x > 0 and math.isfinite(x) break if not ok: lo, hi = 0.0, 1.0 while _pt(hi, df) < p and hi < 1e8: hi *= 2 for _ in range(200): if hi-lo <= 1e-12*max(1, hi): break mid = (lo+hi)/2 if _pt(mid, df) < p: lo = mid else: hi = mid x = (lo+hi)/2 return x def _pnt(t, df, delta): """Noncentral t CDF, algorithm AS 243 (Lenth 1989), following R's pnt().""" if delta == 0: return _pt(t, df) tt, dl, negdel = t, delta, False if t < 0: if delta > 40: return 0.0 negdel, tt, dl = True, -t, -delta if df > 4e5 or dl*dl > 1415.4: s = 1/(4*df) lower = _cdf((tt*(1-s)-dl)/math.sqrt(1+tt*tt*2*s)) return 1-lower if negdel else lower x = tt*tt/(tt*tt+df) tnc = 0.0 if x > 0: lam = dl*dl p = .5*math.exp(-.5*lam) if p == 0: return 1.0 if negdel else 0.0 q = math.sqrt(2/math.pi)*p*dl s = .5-p if s < 1e-7: s = -.5*math.expm1(-.5*lam) a, b = .5, .5*df rxb = (1-x)**b albeta = .5*math.log(math.pi)+math.lgamma(b)-math.lgamma(.5+b) xodd, godd = _ibeta(x, a, b), 2*rxb*math.exp(a*math.log(x)-albeta) tnc = b*x xeven = tnc if tnc < sys.float_info.epsilon else 1-rxb geven = tnc*rxb tnc = p*xodd+q*xeven for it in range(1, 1001): a += 1 xodd -= godd xeven -= geven godd *= x*(a+b-1)/a geven *= x*(a+b-.5)/(a+.5) p *= lam/(2*it) q *= lam/(2*it+1) tnc += p*xodd+q*xeven s -= p if s < -1e-10: break if s <= 0 and it > 1: break if abs(2*s*(xodd-godd)) < 1e-12: break tnc = min(tnc+_cdf(-dl), 1.0) return 1-tnc if negdel else tnc def _t_power(d, n): """Two-sided alpha = .05 power, two equal groups of n (df = 2n-2, ncp = d*sqrt(n/2)); n may be fractional.""" df, ncp = 2*n-2, d*math.sqrt(n/2) crit = _qt(.975, df) return 1 - _pnt(crit, df, ncp) + _pnt(-crit, df, ncp) def _required_per_group(d): """Smallest whole n per group reaching 80% power: start at the normal answer, step to the exact one.""" n = max(2, math.ceil(2*(2.8016/d)**2)) while _t_power(d, n) < .8 and n < 1e7: n += 1 while n > 2 and _t_power(d, n-1) >= .8: n -= 1 return n def simulate_shortform_power(k_full=30, k_short=6, a_mean=1.6, target_zone=(1.0, 2.5), d_true=.4, n_total=250, penalty=.72, boost=1.1): values = (k_full, k_short, a_mean, *target_zone, d_true, n_total, penalty, boost) if not all(math.isfinite(x) for x in values): raise ValueError('All parameters must be finite.') if not (isinstance(k_full, int) and 2 <= k_full <= 200 and isinstance(k_short, int) and 1 <= k_short <= k_full and isinstance(n_total, int) and 4 <= n_total <= 1000000 and .2 <= a_mean <= 4 and -3 <= target_zone[0] < target_zone[1] <= 3 and 0 < d_true <= 3 and 0 < penalty <= 1 and .5 <= boost <= 1.5): raise ValueError('Invalid scenario parameters.') def seq(lo, hi, n): return [(lo + hi) / 2] if n == 1 else [lo + (hi-lo)*i/(n-1) for i in range(n)] def information(theta, a, center): # Cumulative category probabilities and their derivatives give item information. p = [1] + [1/(1+math.exp(-a*(theta-center-t))) for t in (-.9,-.3,.3,.9)] + [0] dp = [a*x*(1-x) for x in p] return sum((dp[i]-dp[i+1])**2/max(p[i]-p[i+1], 1e-12) for i in range(5)) theta = seq(-3, 3, 61) idx = [i for i,t in enumerate(theta) if target_zone[0]-1e-9 <= t <= target_zone[1]+1e-9] if not idx: raise ValueError('Target zone must contain a 0.1-spaced grid point.') specs = {'full': (seq(-2.2,2.2,k_full), a_mean, 1), 'naive': (seq(-.4,.4,k_short), a_mean, penalty), 'targeted': (seq(target_zone[0]-.3,target_zone[1]-.1,k_short), a_mean*boost, 1)} forms = {} for key, (centers,a,weight) in specs.items(): info = [weight*sum(information(t,a,b) for b in centers) for t in theta] # Precision proxy, not empirical or population reliability. rho = sum(info[i]/(info[i]+1) for i in idx)/len(idx) observed = d_true*math.sqrt(rho) forms[key] = dict(info=info, reliability=rho, d=observed, sem=sum(1/math.sqrt(info[i]) for i in idx)/len(idx), power=_t_power(observed, n_total/2), requiredN=2*_required_per_group(observed)) return dict(theta=theta, forms=forms, burden=100*(1-k_short/k_full), extra=forms['naive']['requiredN']-forms['targeted']['requiredN']) if __name__ == '__main__': print(json.dumps(simulate_shortform_power(), indent=2))