brier_index_properness.py

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from scipy.optimize import minimize_scalar
from scipy.stats import binom


def expected_brier_index(f: float, p: float, n: int) -> float:
    """Compute E[Brier Index] for a forecaster who submits forecast f
    on each of n i.i.d. binary questions with true probability p.

    The expectation is taken over the binomial distribution of the
    number of questions that resolve to 1."""
    value = 0
    for k in range(0, n + 1):
        prob_k = binom.pmf(k, n, p)
        avg_brier_score = k / n * (1 - f) ** 2 + (1 - k / n) * (0 - f) ** 2
        value += prob_k * (1 - np.sqrt(avg_brier_score)) * 100
    return value


def find_optimal_forecast(p: float, n: int) -> dict:
    """Find the forecast f that maximizes E[Brier Index] for given p and n."""
    res = minimize_scalar(
        lambda f: -expected_brier_index(f, p, n), bounds=(0, 1), method="bounded"
    )
    f_star = res.x
    expected_brier_index_star = -res.fun
    expected_brier_index_truthful = expected_brier_index(p, p, n)
    return {
        "f_star": f_star,
        "expected_brier_index_star": expected_brier_index_star,
        "expected_brier_index_truthful": expected_brier_index_truthful,
    }


def main():
    p = 0.7  # True probability
    n_values = [1, 10, 30]  # Number of questions
    f_grid = np.linspace(0, 1, 100)  # Grid of forecasts to evaluate

    # Calculate the optimal forecast and expected Brier Index for each n
    optimal_forecasts = []
    for n in n_values:
        optimal_forecasts.append({"n": n, **find_optimal_forecast(p, n)})
    df_optimal = pd.DataFrame(optimal_forecasts)
    df_optimal.to_csv("./optimal_forecasts.csv", index=False)

    # Plot the expected Brier Index as a function of the forecast for each n

    # Evalute the objective function for a grid of forecasts and each n
    res = []
    for f in f_grid:
        for n in n_values:
            res.append(
                {"f": f, "expected_brier_index": expected_brier_index(f, p, n), "n": n}
            )
    df = pd.DataFrame(res)
    df.to_csv("expected_brier_index.csv", index=False)

    # Plot the objective function for different n
    sns.set_theme()
    plt.figure(figsize=(10, 6))
    for n in n_values:
        mask = df["n"] == n
        df_temp = df[mask].copy()
        plt.plot(df_temp["f"], df_temp["expected_brier_index"], label=f"n={n}")
        # Calculate the optimal forecast & annotate
        mask = df_optimal["n"] == n
        opt = df_optimal[mask].iloc[0]
        plt.scatter(opt["f_star"], opt["expected_brier_index_star"])
    plt.axvline(p, color="gray", linestyle="--", label="Truth (p)")
    plt.xlabel("Forecast $f$")
    plt.ylabel("Expected Brier Index")
    plt.title(f"Expected Brier Index vs Forecast $f$ (p={p})")
    plt.legend()
    plt.savefig("./expected_brier_index_plot.png", dpi=300, bbox_inches="tight")
    plt.show()


if __name__ == "__main__":
    main()