1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6EPS = 1e-3
  7
  8
  9def bts_score(p_hat, y, eps=EPS):
 10    return math.log((p_hat + eps) / (y + eps))
 11
 12
 13def prediction_frequency_scaling():
 14    # At fixed y, changing p changes the score by exactly
 15    # log((p1+eps)/(p0+eps)); without epsilon this is log(p1/p0).
 16    y, p0 = 0.20, 0.25
 17    rows = []
 18    for rho in [0.5, 0.8, 1.0, 1.5, 2.0]:
 19        p1 = min(0.95, p0 * rho)
 20        observed = bts_score(p1, y) - bts_score(p0, y)
 21        predicted = math.log((p1 + EPS) / (p0 + EPS))
 22        rows.append(dict(rho=rho, observed=observed, predicted=predicted,
 23                         abs_error=abs(observed - predicted)))
 24    return rows
 25
 26
 27def monotonicity_sweep():
 28    # The derivative with respect to p_hat is 1/(p_hat+eps)>0.
 29    y = 0.35
 30    ps = np.linspace(0.02, 0.98, 25)
 31    rs = np.array([bts_score(float(p), y) for p in ps])
 32    return dict(min_adjacent_delta=float(np.min(np.diff(rs))),
 33                p_at_max=float(ps[np.argmax(rs)]), r_at_min=float(rs[0]),
 34                r_at_max=float(rs[-1]))
 35
 36
 37def variance_scaling(rng):
 38    # For K~Binomial(G,p), delta method predicts
 39    # Var[log((K/G+eps)/(y+eps))] ~= p(1-p)/(G*(p+eps)^2).
 40    p, y, trials = 0.30, 0.20, 50000
 41    rows = []
 42    for G in [8, 16, 32, 64, 128, 256]:
 43        k = rng.binomial(G, p, size=trials)
 44        scores = np.log((k / G + EPS) / (y + EPS))
 45        observed = float(np.var(scores, ddof=1))
 46        predicted = p * (1-p) / (G * (p + EPS)**2)
 47        rows.append(dict(G=G, observed=observed, predicted=predicted,
 48                         observed_times_G=observed*G,
 49                         predicted_times_G=predicted*G,
 50                         ratio=observed/predicted))
 51    return rows
 52
 53
 54def group_reward(answer, predicted, answers, eps=EPS):
 55    counts = np.bincount(answers, minlength=2) / len(answers)
 56    return np.log((counts[answer] + eps) / (predicted + eps))
 57
 58
 59def policy_simulation(rng, rounds=4000, G=16):
 60    # Two answer policies: honest follows latent truth; sycophantic follows
 61    # the user pressure. Each completion predicts the group's answer rate.
 62    # BTS is applied as a group-relative reward; baseline is majority/agreement.
 63    # This is a mechanism simulation, not neural-network training.
 64    results = {}
 65    for method in ["bts", "agreement"]:
 66        for condition in ["neutral", "pressure"]:
 67            correct = flips = 0
 68            reward_var = []
 69            for _ in range(rounds):
 70                truth = int(rng.random() < 0.5)
 71                pressure = truth ^ 1 if condition == "pressure" else truth
 72                # A mixed population: honest samples have 0.78 truth fidelity,
 73                # sycophantic samples have 0.84 pressure fidelity.
 74                honest = rng.random(G) < 0.55
 75                answers = np.where(honest,
 76                                  truth ^ (rng.random(G) > 0.78),
 77                                  pressure ^ (rng.random(G) > 0.84)).astype(int)
 78                # Prediction reports are calibrated to each respondent's rule.
 79                pred = np.where(honest, 0.55 + 0.23 * truth,
 80                                0.55 + 0.29 * pressure).astype(float)
 81                # pred is probability of answer 1; convert for each answer.
 82                ys = np.where(answers == 1, pred, 1-pred)
 83                if method == "bts":
 84                    scores = np.array([group_reward(int(a), float(y), answers)
 85                                       for a, y in zip(answers, ys)])
 86                    reward_var.append(float(np.var(scores)))
 87                    # Reward-weighted soft selection of a policy for the next
 88                    # prompt: BTS favors reports surprising relative to their
 89                    # own prediction, rather than raw consensus.
 90                    w = np.exp(np.clip(scores - scores.max(), -10, 10))
 91                    chosen = int(answers[rng.choice(G, p=w/w.sum())])
 92                else:
 93                    reward = np.bincount(answers, minlength=2) / G
 94                    reward_var.append(float(np.var(reward[answers])))
 95                    chosen = int(np.argmax(np.bincount(answers, minlength=2)))
 96                correct += chosen == truth
 97                flips += condition == "pressure" and chosen != truth
 98            results[f"{method}_{condition}"] = dict(accuracy=correct/rounds,
 99                                                       pressure_flip_rate=flips/rounds,
100                                                       mean_reward_variance=float(np.mean(reward_var)))
101    return results
102
103
104def main():
105    rng = np.random.default_rng(2735)
106    out = {
107        "prediction_1_log_scaling": prediction_frequency_scaling(),
108        "prediction_2_monotonicity": monotonicity_sweep(),
109        "prediction_3_variance_1_over_G": variance_scaling(rng),
110        "policy_simulation": policy_simulation(rng),
111    }
112    Path("results.json").write_text(json.dumps(out, indent=2))
113    print(json.dumps(out, indent=2))
114
115
116if __name__ == "__main__":
117    main()