#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ ============================================================================= Replication code for "Criterion Uncertainty and the Limits of Framing: Identification and Parameter-Free Tests of a Quantum-Probabilistic Choice Model" Chansu Lim ============================================================================= The script reproduces every number, table and figure reported in Section 8 of the manuscript. It is written in *pure Python 3* and uses ONLY the standard library (math, random, statistics). No third-party package is required, so the results can be reproduced on any machine with a stock CPython >= 3.12 interpreter (random.binomialvariate is used). Usage: python qcct_simulation.py Outputs (written to the working directory): qcct_generated.tex LaTeX macros holding all table bodies and all figure coordinates used by the manuscript qcct_results.txt plain-text log of every reported number Fixed seed: 20260619 ============================================================================= """ import math import random import statistics import sys # --------------------------------------------------------------------------- # 0. Global configuration # --------------------------------------------------------------------------- SEED = 20260619 R_REPS_OPT = 5000 # replications per design (state recovery, specification test) R_REPS_FAST = 5000 # replications per design (closed-form test statistics) SAMPLE_SIZES = (200, 500, 2000) # Criterion-conditional logit kernels: probability of alternative 1 given # criterion m. L1 = logistic(1.4), L2 = logistic(-0.8). def logistic(z): return 1.0 / (1.0 + math.exp(-z)) L1 = logistic(1.4) # 0.802184... L2 = logistic(-0.8) # 0.310026... CGAP = L1 - L2 # kernel separation c = L1 - L2 # Data-generating criterion state (pure): Bloch vector GAMMA_TRUE = 1.25 DELTA_TRUE = 0.35 # Coplanar (real) context rotations used in Study A PHIS = (-1.20, -0.60, 0.00, 0.60, 1.20) # Terminal choice effect used in Study C (diagonal in the reference basis) E_HI, E_LO = 0.82, 0.27 # One-sided / two-sided normal critical values and chi-square 95% quantiles Z95_ONE = 1.6448536269514722 # NormalDist().inv_cdf(0.95) Z95_TWO = 1.959963984540054 # NormalDist().inv_cdf(0.975) CHI2_95 = {1: 3.841458820694124, 2: 5.991464547107979, 3: 7.814727903251179, 4: 9.487729036781154} EPS = 1e-12 def clip(x, lo, hi): return lo if x < lo else (hi if x > hi else x) # --------------------------------------------------------------------------- # 1. Qubit primitives (Bloch representation) # --------------------------------------------------------------------------- def bloch_of_pure(gamma, delta): """Bloch vector of cos(gamma/2)|1> + e^{i delta} sin(gamma/2)|2>.""" return (math.sin(gamma) * math.cos(delta), math.sin(gamma) * math.sin(delta), math.cos(gamma)) def dot3(u, v): return u[0] * v[0] + u[1] * v[1] + u[2] * v[2] def cross3(u, v): return (u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]) def q_criterion(r, n): """Born probability of criterion 1 for state r under context direction n.""" return 0.5 * (1.0 + dot3(r, n)) def choice_prob_from_q(q): """P(j = 1 | context) = q L1 + (1 - q) L2.""" return L2 + CGAP * q # --- exact 2x2 complex-matrix machinery, used only to validate the closed # forms derived in the paper (Theorems 12 and 15). ----------------------- def rho_from_bloch(r): """density matrix (I + r . sigma)/2 as a 2x2 tuple of complex numbers.""" x, y, z = r return ((0.5 * (1 + z), 0.5 * complex(x, -y)), (0.5 * complex(x, y), 0.5 * (1 - z))) def projectors(n): """Rank-one projectors of the +-1 observable with Bloch direction n.""" x, y, z = n p_plus = ((0.5 * (1 + z), 0.5 * complex(x, -y)), (0.5 * complex(x, y), 0.5 * (1 - z))) p_minus = ((0.5 * (1 - z), -0.5 * complex(x, -y)), (-0.5 * complex(x, y), 0.5 * (1 + z))) return p_plus, p_minus def mat_mul(a, b): return ((a[0][0] * b[0][0] + a[0][1] * b[1][0], a[0][0] * b[0][1] + a[0][1] * b[1][1]), (a[1][0] * b[0][0] + a[1][1] * b[1][0], a[1][0] * b[0][1] + a[1][1] * b[1][1])) def mat_add(a, b): return ((a[0][0] + b[0][0], a[0][1] + b[0][1]), (a[1][0] + b[1][0], a[1][1] + b[1][1])) def trace(a): return a[0][0] + a[1][1] def dephase(rho, n): """Non-selective state update D_n(rho) = sum_a P_a rho P_a.""" p_plus, p_minus = projectors(n) return mat_add(mat_mul(mat_mul(p_plus, rho), p_plus), mat_mul(mat_mul(p_minus, rho), p_minus)) def terminal_prob(rho): """Tr[E_1 rho] for E_1 = diag(E_HI, E_LO).""" return (E_HI * rho[0][0] + E_LO * rho[1][1]).real # --------------------------------------------------------------------------- # 2. Study A -- recovery of the criterion state and specification test # --------------------------------------------------------------------------- # # Contexts are coplanar rotations in the x-z plane, so only (r_x, r_z) is # identified; it is the natural parameter of Study A. CONTEXT_DIRS_A = tuple((math.sin(p), 0.0, math.cos(p)) for p in PHIS) # With coplanar contexts only the pair (r_x, r_z) enters the likelihood, and # the frame probabilities are AFFINE in it: # p_f(r) = L2 + c (1 + r_x sin phi_f + r_z cos phi_f)/2. # The binomial log-likelihood is therefore globally concave in (r_x, r_z) and # the maximum is unique whenever the design matrix below has rank two. GRAD_A = tuple((0.5 * CGAP * n[0], 0.5 * CGAP * n[2]) for n in CONTEXT_DIRS_A) R_TRUE = bloch_of_pure(GAMMA_TRUE, DELTA_TRUE) RX_TRUE, RZ_TRUE = R_TRUE[0], R_TRUE[2] def probs_A(r_x, r_z): return [choice_prob_from_q(0.5 * (1.0 + r_x * n[0] + r_z * n[2])) for n in CONTEXT_DIRS_A] def loglik_A(r_x, r_z, counts, n_trials): ll = 0.0 for y, p in zip(counts, probs_A(r_x, r_z)): p = clip(p, EPS, 1 - EPS) ll += y * math.log(p) + (n_trials - y) * math.log1p(-p) return ll def loglik_saturated(counts, n_trials): """Unrestricted frame-specific choice probabilities (5 free parameters).""" ll = 0.0 for y in counts: p = clip(y / n_trials, EPS, 1 - EPS) ll += y * math.log(p) + (n_trials - y) * math.log1p(-p) return ll def fit_A(counts, n_trials, tol=1e-10, max_iter=60): """Newton-Raphson maximum likelihood in (r_x, r_z) with analytic score and Hessian. The objective is strictly concave, so the iteration converges to the unique maximiser from any interior start. Returns the estimate, the log-likelihood and the observed information matrix.""" r_x, r_z = 0.0, 0.0 for _ in range(max_iter): gx = gz = hxx = hxz = hzz = 0.0 for y, (dx, dz), p in zip(counts, GRAD_A, probs_A(r_x, r_z)): p = clip(p, 1e-9, 1 - 1e-9) s = y / p - (n_trials - y) / (1.0 - p) w = y / p ** 2 + (n_trials - y) / (1.0 - p) ** 2 gx += s * dx gz += s * dz hxx -= w * dx * dx hxz -= w * dx * dz hzz -= w * dz * dz det = hxx * hzz - hxz * hxz if det <= 0: break step_x = -(hzz * gx - hxz * gz) / det step_z = -(hxx * gz - hxz * gx) / det # damped step: never leave the region where all probabilities are valid t = 1.0 for _ in range(40): cand_x, cand_z = r_x + t * step_x, r_z + t * step_z if all(1e-9 < p < 1 - 1e-9 for p in probs_A(cand_x, cand_z)): break t *= 0.5 r_x, r_z = r_x + t * step_x, r_z + t * step_z if abs(t * step_x) + abs(t * step_z) < tol: break # observed information at the optimum hxx = hxz = hzz = 0.0 for y, (dx, dz), p in zip(counts, GRAD_A, probs_A(r_x, r_z)): p = clip(p, 1e-9, 1 - 1e-9) w = y / p ** 2 + (n_trials - y) / (1.0 - p) ** 2 hxx += w * dx * dx hxz += w * dx * dz hzz += w * dz * dz return (r_x, r_z, loglik_A(r_x, r_z, counts, n_trials), (hxx, hxz, hzz)) def study_A(): true_probs = probs_A(RX_TRUE, RZ_TRUE) rows_recovery, rows_size = [], [] for n_trials in SAMPLE_SIZES: rx, rz, gam = [], [], [] cov_x = cov_z = lr_rej = outside = 0 for _ in range(R_REPS_OPT): counts = [random.binomialvariate(n_trials, p) for p in true_probs] r_x, r_z, ll, (ixx, ixz, izz) = fit_A(counts, n_trials) rx.append(r_x) rz.append(r_z) gam.append(math.acos(clip(r_z, -1 + 1e-12, 1 - 1e-12))) if r_x * r_x + r_z * r_z > 1.0: outside += 1 det = ixx * izz - ixz * ixz if det > 0: se_x = math.sqrt(izz / det) se_z = math.sqrt(ixx / det) if abs(r_x - RX_TRUE) <= Z95_TWO * se_x: cov_x += 1 if abs(r_z - RZ_TRUE) <= Z95_TWO * se_z: cov_z += 1 if 2.0 * (loglik_saturated(counts, n_trials) - ll) > CHI2_95[3]: lr_rej += 1 rows_recovery.append(( n_trials, statistics.fmean(rx) - RX_TRUE, math.sqrt(statistics.fmean([(v - RX_TRUE) ** 2 for v in rx])), cov_x / R_REPS_OPT, statistics.fmean(rz) - RZ_TRUE, math.sqrt(statistics.fmean([(v - RZ_TRUE) ** 2 for v in rz])), cov_z / R_REPS_OPT, math.sqrt(statistics.fmean([(v - GAMMA_TRUE) ** 2 for v in gam])), outside / R_REPS_OPT, )) rows_size.append((n_trials, lr_rej / R_REPS_OPT)) # ---- power of the specification test against non-representable frames -- rows_power = [] for eps in (0.02, 0.04, 0.06): row = [eps] for n_trials in SAMPLE_SIZES: alt = [clip(p + eps * s, 0.01, 0.99) for p, s in zip(true_probs, (1, -1, 1, -1, 1))] rej = 0 for _ in range(R_REPS_OPT): counts = [random.binomialvariate(n_trials, p) for p in alt] _, _, ll, _ = fit_A(counts, n_trials) if 2.0 * (loglik_saturated(counts, n_trials) - ll) > CHI2_95[3]: rej += 1 row.append(rej / R_REPS_OPT) rows_power.append(tuple(row)) return true_probs, rows_recovery, rows_size, rows_power # --------------------------------------------------------------------------- # 3. Study B -- the parameter-free purity restriction # --------------------------------------------------------------------------- # # Three mutually unbiased contexts (Bloch directions x, y, z). For any # criterion state, S = sum_f (2 q^f - 1)^2 = |r|^2 <= 1. An unrestricted # frame-specific latent-class model imposes no such bound. MUB_DIRS = ((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)) def purity_statistic(counts, n_trials): """Bias-corrected S-hat and its standard error (delta method).""" s_hat, var_sum = 0.0, 0.0 for y in counts: p_hat = clip(y / n_trials, EPS, 1 - EPS) u_hat = (2.0 * p_hat - L1 - L2) / CGAP v_hat = 4.0 * p_hat * (1.0 - p_hat) / (n_trials * CGAP ** 2) s_hat += u_hat ** 2 - v_hat # bias correction var_sum += 4.0 * u_hat ** 2 * v_hat # leading-order variance return s_hat, math.sqrt(var_sum) if var_sum > 0 else float('inf') def study_B(): designs = [] # H0: criterion states of decreasing purity (S = |r|^2) for label, scale in (("criterion state, pure", 1.0), ("criterion state, mixed", 0.8), ("criterion state, mixed", 0.6)): r = tuple(scale / math.sqrt(3.0) for _ in range(3)) probs = [choice_prob_from_q(q_criterion(r, n)) for n in MUB_DIRS] designs.append((label, scale ** 2, probs, True)) # H1: frame-specific classical weights that no criterion state can produce for s_target in (1.20, 1.50, 2.00): t = math.sqrt(s_target / 3.0) # 2q - 1 = t in every frame q = 0.5 * (1.0 + t) probs = [choice_prob_from_q(q)] * 3 designs.append(("saturated classical frames", s_target, probs, False)) rows = [] for label, s_true, probs, is_null in designs: row = [label, s_true, is_null] for n_trials in SAMPLE_SIZES: rej = 0 for _ in range(R_REPS_FAST): counts = [random.binomialvariate(n_trials, p) for p in probs] s_hat, se = purity_statistic(counts, n_trials) if (s_hat - 1.0) / se > Z95_ONE: rej += 1 row.append(rej / R_REPS_FAST) rows.append(tuple(row)) return rows # --------------------------------------------------------------------------- # 4. Study C -- sequential contexts and the sin(2 eta) order-effect law # --------------------------------------------------------------------------- def order_pair(eta, r): """Exact P(A then B) and P(B then A) for the terminal choice effect.""" a = (0.0, 0.0, 1.0) b = (math.sin(eta), 0.0, math.cos(eta)) rho = rho_from_bloch(r) p_ab = terminal_prob(dephase(dephase(rho, a), b)) p_ba = terminal_prob(dephase(dephase(rho, b), a)) return p_ab, p_ba def order_gap_closed_form(eta, r): """Proposition 5: (1/2) sin(2 eta) d . (n x r).""" a = (0.0, 0.0, 1.0) b = (math.sin(eta), 0.0, math.cos(eta)) axb = cross3(a, b) norm = math.sqrt(dot3(axb, axb)) if norm < 1e-15: return 0.0 n_hat = tuple(c / norm for c in axb) d = (0.0, 0.0, 0.5 * (E_HI - E_LO)) return 0.5 * math.sin(2.0 * eta) * dot3(d, cross3(n_hat, r)) def study_C(): r = bloch_of_pure(GAMMA_TRUE, DELTA_TRUE) # (i) closed form vs. exact operator computation max_dev = 0.0 curve = [] for k in range(0, 121): eta = math.pi / 2.0 * k / 120.0 p_ab, p_ba = order_pair(eta, r) gap_num = p_ab - p_ba gap_cf = order_gap_closed_form(eta, r) max_dev = max(max_dev, abs(gap_num - gap_cf)) curve.append((eta, p_ab, p_ba, gap_num)) # (ii) detection power of a two-sample test for the order effect rows = [] for k, name in ((0.0, r"$0$"), (0.125, r"$\pi/8$"), (0.25, r"$\pi/4$"), (0.375, r"$3\pi/8$"), (0.5, r"$\pi/2$")): eta = math.pi * k p_ab, p_ba = order_pair(eta, r) row = [name, eta, p_ab - p_ba] for n_trials in SAMPLE_SIZES: rej = 0 for _ in range(R_REPS_FAST): y1 = random.binomialvariate(n_trials, p_ab) y2 = random.binomialvariate(n_trials, p_ba) p1, p2 = y1 / n_trials, y2 / n_trials pool = (y1 + y2) / (2.0 * n_trials) se = math.sqrt(max(pool * (1 - pool) * 2.0 / n_trials, EPS)) if abs(p1 - p2) / se > Z95_TWO: rej += 1 row.append(rej / R_REPS_FAST) rows.append(tuple(row)) return max_dev, curve, rows # --------------------------------------------------------------------------- # 5. Study D -- numerical check of the criterion QQ equality (Proposition 6) # --------------------------------------------------------------------------- def study_D(): """Max violation of p_AB(1,2) + p_AB(2,1) = p_BA(1,2) + p_BA(2,1) over a grid of states and context pairs.""" worst = 0.0 for gk in range(1, 8): for dk in range(0, 8): r_pure = bloch_of_pure(math.pi * gk / 8.0, 2 * math.pi * dk / 8.0) for scale in (1.0, 0.7, 0.3): r = tuple(scale * c for c in r_pure) rho = rho_from_bloch(r) for ek in range(1, 8): eta = math.pi * ek / 8.0 a = (0.0, 0.0, 1.0) b = (math.sin(eta), 0.0, math.cos(eta)) pa = projectors(a) pb = projectors(b) lhs = rhs = 0.0 for i in range(2): for j in range(2): if i == j: continue lhs += trace(mat_mul(mat_mul(mat_mul(pa[i], rho), pa[i]), pb[j])).real rhs += trace(mat_mul(mat_mul(mat_mul(pb[i], rho), pb[i]), pa[j])).real worst = max(worst, abs(lhs - rhs)) return worst # --------------------------------------------------------------------------- # 5b. Study E -- out-of-sample prediction and the attainability (entropy) floor # --------------------------------------------------------------------------- # # Five contexts whose Bloch directions span R^3, so the full state is # identified. One context is held out; the state is estimated from the other # four and used to predict the held-out frame. A saturated frame-specific # mixture makes no prediction at all; the context-invariant mixture is the # natural classical benchmark that does. SQ3 = 1.0 / math.sqrt(3.0) CONTEXT_DIRS_E = ((0.0, 0.0, 1.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (SQ3, SQ3, SQ3), (SQ3, -SQ3, SQ3)) # A mixed criterion state: purity is what makes the attainability floor bite. R_SCALE_E = 0.72 R_TRUE_E = tuple(R_SCALE_E * c for c in bloch_of_pure(GAMMA_TRUE, DELTA_TRUE)) LAM1_TRUE = 0.5 * (1.0 + R_SCALE_E) # top eigenvalue of rho def shannon(p): """Binary Shannon entropy in bits.""" p = clip(p, EPS, 1 - EPS) return -(p * math.log2(p) + (1 - p) * math.log2(1 - p)) def probs_E(r, dirs): return [choice_prob_from_q(q_criterion(r, n)) for n in dirs] def det3(m): return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])) def solve3(A, b): """Solve A x = b for a symmetric positive-definite 3x3 A by Cramer's rule.""" d = det3(A) if abs(d) < 1e-14: return None out = [] for k in range(3): M = [[b[i] if j == k else A[i][j] for j in range(3)] for i in range(3)] out.append(det3(M) / d) return tuple(out) def fit_E(counts, n_trials, dirs, max_iter=80, tol=1e-11): """ML estimate of the Bloch vector r in R^3. The frame probabilities are affine in r, so the binomial log-likelihood is concave and Newton-Raphson converges to the unique maximiser whenever the directions span R^3.""" grads = [tuple(0.5 * CGAP * c for c in n) for n in dirs] r = [0.0, 0.0, 0.0] for _ in range(max_iter): g = [0.0, 0.0, 0.0] H = [[0.0] * 3 for _ in range(3)] for y, d, p in zip(counts, grads, probs_E(tuple(r), dirs)): p = clip(p, 1e-9, 1 - 1e-9) s = y / p - (n_trials - y) / (1.0 - p) w = y / p ** 2 + (n_trials - y) / (1.0 - p) ** 2 for a in range(3): g[a] += s * d[a] for b in range(3): H[a][b] += w * d[a] * d[b] # this is minus the Hessian step = solve3(H, g) if step is None: break t = 1.0 for _ in range(40): cand = tuple(r[a] + t * step[a] for a in range(3)) if all(1e-9 < p < 1 - 1e-9 for p in probs_E(cand, dirs)): break t *= 0.5 move = 0.0 for a in range(3): r[a] += t * step[a] move += abs(t * step[a]) if move < tol: break return tuple(r) def study_E(): true_probs = probs_E(R_TRUE_E, CONTEXT_DIRS_E) n_ctx = len(CONTEXT_DIRS_E) rows = [] for n_trials in SAMPLE_SIZES: se_model = se_inv = 0.0 n_fold = 0 lam_hat, floor_ok = [], 0 for _ in range(R_REPS_OPT): counts = [random.binomialvariate(n_trials, p) for p in true_probs] for h in range(n_ctx): idx = [i for i in range(n_ctx) if i != h] sub_dirs = tuple(CONTEXT_DIRS_E[i] for i in idx) sub_counts = [counts[i] for i in idx] r_hat = fit_E(sub_counts, n_trials, sub_dirs) p_model = probs_E(r_hat, (CONTEXT_DIRS_E[h],))[0] p_inv = sum(sub_counts) / (n_trials * len(idx)) se_model += (p_model - true_probs[h]) ** 2 se_inv += (p_inv - true_probs[h]) ** 2 n_fold += 1 # attainability floor implied by the full five-context fit r_full = fit_E(counts, n_trials, CONTEXT_DIRS_E) norm = math.sqrt(sum(c * c for c in r_full)) l1 = 0.5 * (1.0 + min(norm, 1.0)) lam_hat.append(l1) q_obs = [clip((y / n_trials - L2) / CGAP, 0.0, 1.0) for y in counts] q_max = max(max(q, 1.0 - q) for q in q_obs) if l1 + 1e-12 >= q_max: floor_ok += 1 rows.append(( n_trials, math.sqrt(se_model / n_fold), math.sqrt(se_inv / n_fold), statistics.fmean(lam_hat) - LAM1_TRUE, math.sqrt(statistics.fmean([(v - LAM1_TRUE) ** 2 for v in lam_hat])), floor_ok / R_REPS_OPT, )) return true_probs, rows # --------------------------------------------------------------------------- # 6. Reporting # --------------------------------------------------------------------------- def fmt(x, d=3): return ("%." + str(d) + "f") % x def texpow(x): """Format a small positive number as LaTeX scientific notation.""" mant, exp = ("%.1e" % x).split("e") return "%s\\times 10^{%d}" % (mant, int(exp)) def main(): random.seed(SEED) log = [] def say(line=""): print(line) log.append(line) say("=" * 74) say("Criterion uncertainty in discrete choice -- replication output") say("seed = %d ; python = %s" % (SEED, sys.version.split()[0])) say("L1 = %.6f , L2 = %.6f , kernel gap c = %.6f" % (L1, L2, CGAP)) say("=" * 74) # ---------------- Study A ------------------------------------------------ true_probs, rec, size, power = study_A() say("") say("STUDY A -- criterion-state recovery, %d replications" % R_REPS_OPT) say("exact frame choice probabilities: " + ", ".join(fmt(p, 4) for p in true_probs)) say(" n bias(rx) rmse(rx) cov95 bias(rz) rmse(rz) cov95 " "rmse(gamma) outside disc") for row in rec: say(" %-5d %+8.4f %8.4f %6.3f %+8.4f %8.4f %6.3f %10.4f %8.3f" % row) say(" specification test (LR vs unrestricted frames, df = 3, 5%% level)") for n_trials, r_ in size: say(" n = %-5d empirical size = %.3f" % (n_trials, r_)) say(" power against non-representable frame probabilities") for row in power: say(" eps = %.2f : " % row[0] + " ".join("n=%d: %.3f" % (n, v) for n, v in zip(SAMPLE_SIZES, row[1:]))) # ---------------- Study B ------------------------------------------------ rows_B = study_B() say("") say("STUDY B -- purity restriction S <= 1, %d replications" % R_REPS_FAST) say(" design S_true n=200 n=500 n=2000") for row in rows_B: say(" %-28s %6.2f %6.3f %6.3f %6.3f" % (row[0].replace("$", "").replace("\\", ""), row[1], row[3], row[4], row[5])) # ---------------- Study C ------------------------------------------------ max_dev, curve, rows_C = study_C() say("") say("STUDY C -- sequential contexts, %d replications" % R_REPS_FAST) say(" max |exact - closed form| over 121 angles = %.3e" % max_dev) say(" eta gap n=200 n=500 n=2000") for row in rows_C: say(" %-9s %+7.4f %6.3f %6.3f %6.3f" % (row[0].replace("$", "").replace("\\", ""), row[2], row[3], row[4], row[5])) # ---------------- Study D ------------------------------------------------ worst = study_D() say("") say("STUDY D -- criterion QQ equality") say(" max violation over 1176 (state, context-pair) combinations = %.3e" % worst) # ---------------- Study E ------------------------------------------------ probs_E_true, rows_E = study_E() say("") say("STUDY E -- out-of-sample prediction, %d replications" % R_REPS_OPT) say(" five spanning contexts; exact frame probabilities: " + ", ".join(fmt(p, 4) for p in probs_E_true)) say(" true top eigenvalue lambda_1 = %.4f ; entropy floor = %.4f bits" % (LAM1_TRUE, shannon(LAM1_TRUE))) say(" n RMSE(model) RMSE(invariant) bias(lam1) rmse(lam1) floor ok") for row in rows_E: say(" %-5d %10.4f %14.4f %+10.4f %10.4f %8.3f" % row) # ---------------- exact context-response curve for Figure 3 -------------- r_true = bloch_of_pure(GAMMA_TRUE, DELTA_TRUE) resp = [] for k in range(0, 97): phi = -math.pi + 2 * math.pi * k / 96.0 n = (math.sin(phi), 0.0, math.cos(phi)) q_coh = q_criterion(r_true, n) # dephased benchmark: coherence in the reference basis removed q_dep = 0.5 * (1.0 + math.cos(GAMMA_TRUE) * math.cos(phi)) resp.append((phi, choice_prob_from_q(q_coh), choice_prob_from_q(q_dep))) # ---------------- LaTeX artefacts --------------------------------------- with open("qcct_generated.tex", "w", encoding="utf-8") as fh: w = fh.write w("%% Auto-generated by qcct_simulation.py (seed %d). Do not edit.\n" % SEED) w("\n%% ---- Table: parameter recovery -------------------------------\n") w("\\newcommand{\\TabRecovery}{%\n") for row in rec: w("%d & %s & %s & %s & %s & %s & %s & %s \\\\\n" % (row[0], fmt(row[1], 4), fmt(row[2], 4), fmt(row[3]), fmt(row[4], 4), fmt(row[5], 4), fmt(row[6]), fmt(row[7], 4))) w("}\n") w("\n%% ---- Table: specification test size and power ----------------\n") w("\\newcommand{\\TabSpec}{%\n") w("size (data from the model) & %s \\\\\n" % " & ".join(fmt(v) for _, v in size)) for row in power: w("power, $\\varepsilon = %.2f$ & %s \\\\\n" % (row[0], " & ".join(fmt(v) for v in row[1:]))) w("}\n") w("\n%% ---- Table: out-of-sample prediction --------------------------\n") w("\\newcommand{\\TabOOS}{%\n") for row in rows_E: w("%d & %s & %s & %s & %s & %s \\\\\n" % (row[0], fmt(row[1], 4), fmt(row[2], 4), fmt(row[3], 4), fmt(row[4], 4), fmt(row[5]))) w("}\n") w("\n%% ---- Table: purity test --------------------------------------\n") w("\\newcommand{\\TabPurity}{%\n") for row in rows_B: w("%s & %.2f & %s & %s & %s \\\\\n" % (row[0], row[1], fmt(row[3]), fmt(row[4]), fmt(row[5]))) w("}\n") w("\n%% ---- Table: order-effect detection ---------------------------\n") w("\\newcommand{\\TabOrder}{%\n") for row in rows_C: w("%s & %s & %s & %s & %s \\\\\n" % (row[0], fmt(row[2], 4), fmt(row[3]), fmt(row[4]), fmt(row[5]))) w("}\n") w("\n%% ---- Figure coordinates: context response --------------------\n") w("\\newcommand{\\PlotCoherent}{%\n") w(" ".join("(%.4f,%.4f)" % (p[0], p[1]) for p in resp)) w("}\n\\newcommand{\\PlotDephased}{%\n") w(" ".join("(%.4f,%.4f)" % (p[0], p[2]) for p in resp)) w("}\n") w("\n%% ---- Figure coordinates: order effect ------------------------\n") w("\\newcommand{\\PlotOrderAB}{%\n") w(" ".join("(%.4f,%.4f)" % (c[0], c[1]) for c in curve)) w("}\n\\newcommand{\\PlotOrderBA}{%\n") w(" ".join("(%.4f,%.4f)" % (c[0], c[2]) for c in curve)) w("}\n\\newcommand{\\PlotOrderGap}{%\n") w(" ".join("(%.4f,%.4f)" % (c[0], c[3]) for c in curve)) w("}\n") w("\n%% ---- Scalars -------------------------------------------------\n") w("\\newcommand{\\MaxOrderDev}{%s}\n" % texpow(max_dev)) w("\\newcommand{\\MaxQQDev}{%s}\n" % texpow(worst)) w("\\newcommand{\\NRepsOpt}{%d}\n" % R_REPS_OPT) w("\\newcommand{\\NRepsFast}{%s}\n" % f"{R_REPS_FAST:,}") w("\\newcommand{\\KernelOne}{%.3f}\n" % L1) w("\\newcommand{\\KernelTwo}{%.3f}\n" % L2) w("\\newcommand{\\RxTrue}{%.4f}\n" % RX_TRUE) w("\\newcommand{\\RzTrue}{%.4f}\n" % RZ_TRUE) w("\\newcommand{\\RNormTrue}{%.4f}\n" % math.sqrt(RX_TRUE ** 2 + RZ_TRUE ** 2)) w("\\newcommand{\\FrameProbs}{%s}\n" % ", ".join("%.3f" % p for p in true_probs)) w("\\newcommand{\\OutsideDisc}{%.3f}\n" % rec[0][8]) w("\\newcommand{\\LamOneTrue}{%.3f}\n" % LAM1_TRUE) w("\\newcommand{\\EntropyFloor}{%.3f}\n" % shannon(LAM1_TRUE)) w("\\newcommand{\\RScaleE}{%.2f}\n" % R_SCALE_E) w("\\newcommand{\\OrderGapMax}{%.4f}\n" % max(abs(c[3]) for c in curve)) with open("qcct_results.txt", "w", encoding="utf-8") as fh: fh.write("\n".join(log) + "\n") say("") say("wrote qcct_generated.tex and qcct_results.txt") if __name__ == "__main__": main()