Estimator comparison: budget, variance, and accuracy

Business context: Attribution quality must be judged against exact Shapley with realistic player counts and tight call budgets — not from a single seed.

Goal: (1) Benchmark eight players with at most 20 coalition evaluations; (2) sweep budget.fraction up to 0.5 where that cap allows; (3) explain why Monte Carlo can show high error with tiny variance.

Prerequisites: nlp-shap, Python 3.12.

Theory: Approximation methods.

from collections.abc import Callable
from itertools import product

import numpy as np

from nlp_shap import (
    ComplementaryEstimator,
    ExactEstimator,
    MonteCarloEstimator,
    NeymanEstimator,
    ShapleyAggregator,
)
from nlp_shap.domain.coalition import CoalitionMask
from nlp_shap.domain.players import PlayerSet
from nlp_shap.estimation._shared import (
    compute_complementary_num_splits,
    compute_mc_num_samples,
    present_to_mask_int,
)
from nlp_shap.estimation.neyman import _NeymanStep

MAX_EVALUATIONS = 20
SEED_COUNT = 32

Shared benchmark utilities

Technical: Nonlinear characteristic function, exact enumeration reference, and helpers that map a target cap (≤ 20) to the budget.fraction each estimator needs.

PayoffFn = Callable[[tuple[bool, ...]], float]


def interaction_payoff(present: tuple[bool, ...], num_players: int) -> float:
    coefficients = tuple(float(index + 1) for index in range(num_players))
    linear = sum(
        coefficient
        for coefficient, is_present in zip(coefficients, present, strict=True)
        if is_present
    )
    return linear + 0.2 * (linear * linear) / num_players


def exact_reference(num_players: int) -> tuple[list[float], float, tuple[CoalitionMask, ...]]:
    player_set = PlayerSet(player_ids=tuple(f"p{index}" for index in range(num_players)))
    payoff_fn: PayoffFn = lambda present: interaction_payoff(present, num_players)
    all_masks = tuple(
        CoalitionMask.from_sequence(mask)
        for mask in product([False, True], repeat=num_players)
        if not all(mask)
    )
    all_payoffs = [payoff_fn(mask.present) for mask in all_masks]
    exact_values = ExactEstimator().estimate_attributions(
        all_masks,
        all_payoffs,
        ShapleyAggregator(),
    )
    scale = sum(abs(value) for value in exact_values)
    return exact_values, scale, all_masks


def mc_fraction_for_cap(num_players: int, target_cap: int) -> float:
    maximum = (1 << num_players) - 1
    return min(1.0, target_cap / maximum)


def cc_fraction_for_cap(num_players: int, target_cap: int) -> float:
    maximum = (1 << num_players) - 2
    return min(1.0, target_cap / maximum)


def l1_error(estimated: list[float], reference: list[float]) -> float:
    return sum(
        abs(left - right) for left, right in zip(estimated, reference, strict=True)
    )


def run_monte_carlo(
    player_set: PlayerSet,
    exact_values: list[float],
    budget_fraction: float,
    payoff_fn: PayoffFn,
    seed: int,
    include_minimal_masks: bool = True,
) -> float:
    estimator = MonteCarloEstimator()
    masks = list(
        estimator.sample_masks(
            player_set,
            budget_fraction=budget_fraction,
            include_minimal_masks=include_minimal_masks,
            seed=seed,
        )
    )
    payoffs = [payoff_fn(mask.present) for mask in masks]
    values = estimator.estimate_attributions(masks, payoffs, ShapleyAggregator())
    return l1_error(values, exact_values)


def run_complementary(
    player_set: PlayerSet,
    exact_values: list[float],
    budget_fraction: float,
    payoff_fn: PayoffFn,
    seed: int,
) -> float:
    estimator = ComplementaryEstimator()
    masks = list(
        estimator.sample_masks(
            player_set,
            budget_fraction=budget_fraction,
            include_minimal_masks=True,
            seed=seed,
        )
    )
    payoffs = [payoff_fn(mask.present) for mask in masks]
    values = estimator.estimate_attributions(masks, payoffs)
    return l1_error(values, exact_values)


def run_neyman(
    player_set: PlayerSet,
    exact_values: list[float],
    budget_fraction: float,
    payoff_fn: PayoffFn,
    seed: int,
) -> float:
    estimator = NeymanEstimator(use_standard_method=False)
    phase_one = list(
        estimator.sample_masks(
            player_set,
            budget_fraction=budget_fraction,
            include_minimal_masks=False,
            seed=seed,
        )
    )
    payoffs_one = [payoff_fn(mask.present) for mask in phase_one]
    masks = phase_one
    payoffs = payoffs_one
    if estimator._step == _NeymanStep.NEYMAN_ALLOCATION:
        estimator.begin_allocation(phase_one, payoffs_one)
        phase_two = list(estimator.sample_allocation_masks())
        masks = phase_one + phase_two
        payoffs = payoffs_one + [payoff_fn(mask.present) for mask in phase_two]
    values = estimator.estimate_attributions(masks, payoffs)
    return l1_error(values, exact_values)


def evaluate_all(
    num_players: int,
    mc_fraction: float,
    cc_fraction: float,
    seed_count: int,
) -> list[dict[str, object]]:
    player_set = PlayerSet(player_ids=tuple(f"p{index}" for index in range(num_players)))
    exact_values, scale, _ = exact_reference(num_players)
    payoff_fn: PayoffFn = lambda present: interaction_payoff(present, num_players)
    mc_cap = compute_mc_num_samples(num_players, mc_fraction, True)
    cc_cap = compute_complementary_num_splits(num_players, cc_fraction, True)
    rows: list[dict[str, object]] = []
    for label, runner, cap in (
        ("mc", lambda seed: run_monte_carlo(player_set, exact_values, mc_fraction, payoff_fn, seed), mc_cap),
        ("complementary", lambda seed: run_complementary(player_set, exact_values, cc_fraction, payoff_fn, seed), cc_cap),
        ("neyman_cc", lambda seed: run_neyman(player_set, exact_values, cc_fraction, payoff_fn, seed), cc_cap),
    ):
        errors = [runner(seed) for seed in range(seed_count)]
        array = np.asarray(errors, dtype=np.float64)
        rows.append(
            {
                "estimator": label,
                "cap": cap,
                "mean_l1": float(array.mean()),
                "std_l1": float(array.std()),
                "mean_rel": float(array.mean() / scale),
            }
        )
    return rows


def print_rows(title: str, rows: list[dict[str, object]], extra: str = "") -> None:
    print(title)
    if extra:
        print(extra)
    header = f"{'estimator':13s} {'cap':>4s} {'mean L1':>9s} {'std L1':>8s} {'mean rel':>9s}"
    print(header)
    print("-" * len(header))
    for row in rows:
        print(
            f"{row['estimator']:13s} {row['cap']:4d} "
            f"{row['mean_l1']:9.3f} {row['std_l1']:8.3f} {row['mean_rel']:9.3f}"
        )

Eight players, cap ≤ 20 (primary benchmark)

Business: Eight tokens is a realistic sentence slice. We cap coalition evaluations at 20 and ladder toward that maximum.

Technical: At n = 8 a budget.fraction of 0.5 would imply ~127 calls — above the cap — so we convert target caps {8, 12, 16, 20} into equivalent fractions (~0.03–0.08).

NUM_PLAYERS_EIGHT = 8
TARGET_CAPS = (8, 12, 16, 20)

exact_eight, scale_eight, _ = exact_reference(NUM_PLAYERS_EIGHT)
print(
    f"exact Shapley ({NUM_PLAYERS_EIGHT} players): "
    f"{[round(value, 3) for value in exact_eight]}"
)
print(f"reference L1 scale: {scale_eight:.3f}")

for target_cap in TARGET_CAPS:
    assert target_cap <= MAX_EVALUATIONS
    mc_fraction = mc_fraction_for_cap(NUM_PLAYERS_EIGHT, target_cap)
    cc_fraction = cc_fraction_for_cap(NUM_PLAYERS_EIGHT, target_cap)
    rows = evaluate_all(NUM_PLAYERS_EIGHT, mc_fraction, cc_fraction, SEED_COUNT)
    print_rows(
        f"\ntarget cap {target_cap} (mc_fraction={mc_fraction:.3f}, cc_fraction={cc_fraction:.3f})",
        rows,
    )

exact Shapley (8 players): [1.553, 3.112, 4.678, 6.25, 7.828, 9.412, 11.003, 12.6]
reference L1 scale: 56.437

target cap 8 (mc_fraction=0.031, cc_fraction=0.031)
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc               9    51.300    0.000     0.909
complementary   16    54.688    6.080     0.969
neyman_cc       16    35.682   10.753     0.632

target cap 12 (mc_fraction=0.047, cc_fraction=0.047)
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc              12    51.240    0.087     0.908
complementary   16    54.688    6.080     0.969
neyman_cc       16    35.682   10.753     0.632

target cap 16 (mc_fraction=0.063, cc_fraction=0.063)
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc              16    51.083    0.187     0.905
complementary   16    54.688    6.080     0.969
neyman_cc       16    35.682   10.753     0.632

target cap 20 (mc_fraction=0.078, cc_fraction=0.079)
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc              20    50.961    0.213     0.903
complementary   20    48.080   10.109     0.852
neyman_cc       20    36.634    8.686     0.649

Why Monte Carlo shows high error but tiny variance

Business: A narrow error band is not proof of accuracy — MC can be confidently wrong.

Technical:

1. include_minimal_masks=True emits the same empty/singleton prefix on every seed. When the cap is ≤ n + 1, the sample is deterministic.

2. ShapleyAggregator.aggregate_from_marginals skips Shapley terms unless both S and S ∪ {i} were sampled — most terms never enter the sum.

def shapley_marginal_coverage(
    masks: list[CoalitionMask],
    num_players: int,
) -> tuple[int, int]:
    present_lookup = {present_to_mask_int(mask.present) for mask in masks}
    resolved = 0
    possible = 0
    for player_index in range(num_players):
        player_bit = 1 << player_index
        for mask in masks:
            mask_int = present_to_mask_int(mask.present)
            if mask_int & player_bit:
                continue
            possible += 1
            if (mask_int | player_bit) in present_lookup:
                resolved += 1
    return resolved, possible


player_set_eight = PlayerSet(
    player_ids=tuple(f"p{index}" for index in range(NUM_PLAYERS_EIGHT))
)
payoff_eight: PayoffFn = lambda present: interaction_payoff(present, NUM_PLAYERS_EIGHT)
analysis_fraction = mc_fraction_for_cap(NUM_PLAYERS_EIGHT, 12)

prefix_a = [
    mask.present
    for mask in list(
        MonteCarloEstimator().sample_masks(
            player_set_eight, analysis_fraction, True, seed=0
        )
    )[: NUM_PLAYERS_EIGHT + 1]
]
prefix_b = [
    mask.present
    for mask in list(
        MonteCarloEstimator().sample_masks(
            player_set_eight, analysis_fraction, True, seed=99
        )
    )[: NUM_PLAYERS_EIGHT + 1]
]
print(f"minimal prefix identical across seeds: {prefix_a == prefix_b}")

sample_masks = list(
    MonteCarloEstimator().sample_masks(
        player_set_eight,
        budget_fraction=analysis_fraction,
        include_minimal_masks=True,
        seed=0,
    )
)
resolved, possible = shapley_marginal_coverage(sample_masks, NUM_PLAYERS_EIGHT)
print(
    f"Shapley marginals resolved: {resolved}/{possible} "
    f"({100 * resolved / possible:.1f}%) at cap={len(sample_masks)}"
)

for include_minimal_masks in (True, False):
    errors = [
        run_monte_carlo(
            player_set_eight,
            exact_eight,
            analysis_fraction,
            payoff_eight,
            seed,
            include_minimal_masks=include_minimal_masks,
        )
        for seed in range(SEED_COUNT)
    ]
    cap = compute_mc_num_samples(
        NUM_PLAYERS_EIGHT, analysis_fraction, include_minimal_masks
    )
    print(
        f"include_minimal_masks={include_minimal_masks} cap={cap} "
        f"mean L1={np.mean(errors):.3f} std L1={np.std(errors):.3f}"
    )

minimal prefix identical across seeds: True
Shapley marginals resolved: 9/73 (12.3%) at cap=12
include_minimal_masks=True cap=12 mean L1=51.240 std L1=0.087
include_minimal_masks=False cap=12 mean L1=56.358 std L1=0.071

Fraction sweep 0.1–0.5 at cap ≤ 20 (five players)

Business: Product configs often speak in budget.fraction, not coalition counts. At eight players, fraction = 0.5 would exceed 20 calls, so we repeat the sweep on five players where 0.1–0.5 all stay within the cap.

Technical: Same nonlinear game and estimators; exact reference recomputed for n = 5.

NUM_PLAYERS_FIVE = 5
BUDGET_FRACTIONS = (0.1, 0.2, 0.3, 0.4, 0.5)

exact_five, scale_five, _ = exact_reference(NUM_PLAYERS_FIVE)
print(
    f"exact Shapley ({NUM_PLAYERS_FIVE} players): "
    f"{[round(value, 3) for value in exact_five]}"
)

for budget_fraction in BUDGET_FRACTIONS:
    mc_cap = compute_mc_num_samples(NUM_PLAYERS_FIVE, budget_fraction, True)
    cc_cap = compute_complementary_num_splits(NUM_PLAYERS_FIVE, budget_fraction, True)
    assert mc_cap <= MAX_EVALUATIONS
    assert cc_cap <= MAX_EVALUATIONS
    rows = evaluate_all(
        NUM_PLAYERS_FIVE,
        budget_fraction,
        budget_fraction,
        SEED_COUNT,
    )
    print_rows(f"\nbudget.fraction = {budget_fraction:.1f}", rows)

print("\nexact reference: L1 = 0.000 (by definition)")

exact Shapley (5 players): [1.168, 2.352, 3.552, 4.768, 6.0]

budget.fraction = 0.1
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc               6    14.400    0.000     0.807
complementary   10    15.754    3.250     0.883
neyman_cc       10    10.595    2.017     0.594

budget.fraction = 0.2
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc               6    14.400    0.000     0.807
complementary   10    15.754    3.250     0.883
neyman_cc       10    10.595    2.017     0.594

budget.fraction = 0.3
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc               9    13.901    0.316     0.779
complementary   10    15.754    3.250     0.883
neyman_cc       10    10.595    2.017     0.594

budget.fraction = 0.4
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc              12    12.844    0.446     0.720
complementary   12    12.933    3.600     0.725
neyman_cc       12    11.357    1.999     0.637

budget.fraction = 0.5
estimator      cap   mean L1   std L1  mean rel
-----------------------------------------------
mc              15    11.754    0.465     0.659
complementary   14    11.142    2.085     0.625
neyman_cc       14    12.589    1.983     0.706

exact reference: L1 = 0.000 (by definition)

Reading the results

Business: On eight players with ≤ 20 calls, Neyman-CC and complementary track exact Shapley more closely than MC. Do not trust MC’s small variance — it reflects a fixed, incomplete sample, not convergence.

Technical: On the five-player fraction ladder, Neyman leads at low fractions; complementary catches up near 0.5. MC improves with cap but remains biased while aggregate_from_marginals drops unresolved terms.