1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6
  7def sigmoid(x):
  8    return 1.0 / (1.0 + np.exp(-np.clip(x, -50.0, 50.0)))
  9
 10
 11def q_values(p_m, eps, beta, p_a):
 12    q_h = (1 - p_m) * eps + p_m * (1 - eps)
 13    return q_h, (q_h + beta * p_a) / (1 + beta)
 14
 15
 16def update_pm(p_m, rho, kappa, delta_pi):
 17    return (1 - rho) * p_m + rho * sigmoid(kappa * delta_pi)
 18
 19
 20def map_decision(reports, qs):
 21    reports = np.asarray(reports)
 22    qs = np.clip(np.asarray(qs), 1e-3, 1 - 1e-3)
 23    score = np.sum(np.where(reports == 1, np.log((1-qs)/qs),
 24                            np.log(qs/(1-qs))))
 25    return 1 if score >= 0 else -1
 26
 27
 28def implied_pa_threshold(p_m, eps, beta):
 29    """P_a solving q=1/2 under the stated q formula."""
 30    qh, _ = q_values(p_m, eps, beta, 0.0)
 31    return (0.5 * (1 + beta) - qh) / beta
 32
 33
 34def math_checks():
 35    eps, beta = .12, .35
 36    p = 0.0
 37    for _ in range(300):
 38        p = update_pm(p, .2, 3., 1.)
 39    qh0, q0 = q_values(0, eps, beta, .7)
 40    qh1, q1 = q_values(1, eps, beta, .7)
 41    return {
 42        "bounded_dynamics": bool(0 <= p <= 1),
 43        "q_formula_endpoints": bool(abs(qh0-eps) < 1e-12 and abs(qh1-(1-eps)) < 1e-12),
 44        "q_at_pa_half_p0": q_values(0, eps, beta, .5)[1],
 45        "q_at_pa_half_p1": q_values(1, eps, beta, .5)[1],
 46        "implied_pa_threshold_p0": implied_pa_threshold(0, eps, beta),
 47        "implied_pa_threshold_p1": implied_pa_threshold(1, eps, beta),
 48        "claimed_universal_threshold_holds": bool(
 49            abs(q_values(0, eps, beta, .5)[1]-.5) < 1e-10 and
 50            abs(q_values(1, eps, beta, .5)[1]-.5) < 1e-10),
 51    }
 52
 53
 54def one_round(rng, n_honest, beta, eps, pa, pm, method, adaptive):
 55    n_bad = int(round(beta*n_honest))
 56    truth = 1 if rng.random() < .5 else -1
 57    sensed = np.where(rng.random(n_honest) < eps, -truth, truth)
 58    honest = np.where(rng.random(n_honest) < pm, -sensed, sensed)
 59    bad = np.where(rng.random(n_bad) < pa, -truth, truth)
 60    reports = np.r_[honest, bad]
 61    q = q_values(pm, eps, beta, pa)[1]
 62    if method == 'majority':
 63        pred = 1 if reports.sum() >= 0 else -1
 64    elif method == 'map':
 65        pred = map_decision(reports, np.full(reports.size, q))
 66    else:
 67        raise ValueError(method)
 68    observed_error = np.mean(reports != truth)
 69    # Simple bounded-rational conformity proxy, included only to exercise the
 70    # proposed state update; it is not claimed to reproduce the paper's payoff.
 71    next_pm = update_pm(pm, .25, 4., 2*(observed_error-.5)) if adaptive else pm
 72    return int(pred == truth), next_pm, q, observed_error
 73
 74
 75def simulate(rng, pa, beta, eps, method, adaptive, rounds=20):
 76    pm, acc, qs, observed = 0., [], [], []
 77    for _ in range(rounds):
 78        a, pm, q, e = one_round(rng, 20, beta, eps, pa, pm, method, adaptive)
 79        acc.append(a); qs.append(q); observed.append(e)
 80    return np.mean(acc), np.mean(qs), pm, np.mean(observed)
 81
 82
 83def run_experiment(seed=7):
 84    rng = np.random.default_rng(seed)
 85    rows = []
 86    # Two beta values make dependence on attacker population explicit.
 87    for beta in (.25, 1.0, 2.0):
 88        for eps in (.05, .20):
 89            for pa in np.linspace(.1, .9, 9):
 90                for method, adaptive in [('majority', False), ('map', False),
 91                                         ('map', True)]:
 92                    vals = [simulate(rng, float(pa), beta, eps, method, adaptive)
 93                            for _ in range(60)]
 94                    rows.append({
 95                        'beta': beta, 'eps': eps, 'p_a': float(pa),
 96                        'method': method, 'adaptive': adaptive,
 97                        'accuracy': float(np.mean([v[0] for v in vals])),
 98                        'mean_q': float(np.mean([v[1] for v in vals])),
 99                        'final_pm': float(np.mean([v[2] for v in vals])),
100                        'observed_error': float(np.mean([v[3] for v in vals]))})
101    return rows
102
103
104def main():
105    out = {'math_checks': math_checks(), 'results': run_experiment()}
106    Path('results.json').write_text(json.dumps(out, indent=2))
107    print(json.dumps(out['math_checks'], indent=2))
108    for r in out['results']:
109        if r['eps'] == .05 and r['method'] in ('majority', 'map') and not r['adaptive']:
110            print(r)
111
112
113if __name__ == '__main__':
114    main()