# Experiment 1: fresh-sample and historical-mean gradients

We solve the same stochastic problem using two gradient estimators:

\[
\min_\theta\;\mathbb E[\tfrac12(\theta-X)^2+\log\cosh\theta],
\qquad X\sim N(m,1).
\]
The algorithms observe samples, not the population mean. The default experiment uses `m=0`, so `theta*=0` and `H=2`.

- Fresh gradient: `theta + tanh(theta) - X_t`.
- Historical gradient: `theta + tanh(theta) - mean(X_1,...,X_t)`.

Both algorithms receive the same samples within each replication and use the same initial value and step schedule. Replications are independent. The PR average is `(theta_1+...+theta_T)/T`, including pre-update iterates, with no discarded warm-up. Each histogram entry is one independent replication of `sqrt(T)*(PR-theta*)`.

Defaults: 10,000 replications; checkpoints `T=1000,5000,20000,50000,100000,200000`; `eta_t=0.5*t^(-0.65)`; initial iterate 2; seed 20261016. The simulation vectorizes across replications and retains only checkpoint draws.

## Run

From this directory, using the local environment created for the experiment:

```sh
.venv/bin/python pr_histograms.py
```

For a fresh Python environment:

```sh
python3 -m venv .venv
.venv/bin/python -m pip install -r requirements.txt
.venv/bin/python pr_histograms.py
```

Parameters are configurable through `--reps`, `--horizons`, `--seed`, `--alpha`, `--step-scale`, `--initial`, `--mean`, `--bins`, and `--output`.

## Outputs

Saved under `results/experiment_1/`:

- `pr_histograms.png` and `.pdf`: shared-axis histograms at the largest horizon, with theoretical normal densities and empirical moments.
- `variance_convergence.png` and `.pdf`: empirical variances versus horizon, with asymptotic reference lines and approximate Monte Carlo error bars.
- `scaled_pr_draws.npz`: every checkpoint's scaled PR draws, allowing independent verification or replotting.
- `summary.json` and `summary.csv`: configuration, theoretical values, empirical moments, Monte Carlo standard errors, and gradient-oracle MSE at each checkpoint.

At `m=0`, the predicted total PR variances are 1/4 (fresh samples) and 1/2 (historical mean). Their difference compares two different noise structures. The single-step oracle MSEs at a common fixed parameter are 1 and 1/t, respectively. These are compared with errors relative to the population mean, not with the difference between the two algorithm trajectories. Variance error bars quantify Monte Carlo uncertainty and do not account for finite-horizon bias.

## Experiment 2: estimator-gain sweep

Run:

```sh
.venv/bin/python gain_sweep.py
```

This reuses the objective and SGD settings from Experiment 1. The auxiliary estimator starts at zero and uses the exact recursion `mhat_t=(1-a/t)*mhat_{t-1}+(a/t)*X_t`, including early iterations. The gain grid is `0.6,0.75,1,1.25,1.5,2,3,5,8`. There are 5,000 independent replications per gain, paired by sharing each replication's sample stream across gains. Checkpoints are 50,000 and 200,000 iterations. The seed is 20261017. The gain `a=1` is the running-average method of Experiment 1; this experiment uses its own replication seed.

Outputs are in `results/experiment_2/`: `gain_sweep.pdf`, `gain_sweep.png`, `gain_draws.npz`, `summary.json`, and `summary.csv`. The figure compares:

- `T * E[(mhat_T-m)^2]`, with asymptotic curve `a^2/(2a-1)`;
- the empirical variance of `sqrt(T)*(PR-theta*)`, with asymptotic curve `2a/((2a-1)*H^2)` and lower-bound line `1/H^2`.

The first curve has its asymptotic minimum at `a=1`. The second decreases with the gain. Error bars show 1.96 Monte Carlo standard errors. The two displayed horizons expose finite-horizon effects, especially near the threshold `a=1/2`; the variance limits are not assumed to hold exactly at finite iteration counts. A gain above one can produce signed weights during the finite initial interval; the implementation follows the paper's unshifted `a/t` recursion exactly.

## Experiment 3: attaining the downstream variance bound

```sh
.venv/bin/python attainment_histograms.py
```

The theoretical PR variances are printed and saved in `results/experiment_3/theory.json` before simulation: 1/2 for the running average, 4/15 for the 8/t recursion, and 1/4 for the t^(-0.6) recursion. All three estimators start at zero and follow the specified gain from the first step; the polynomial recursion has first gain one and thus exact weight normalization. Other SGD settings match Experiment 1. Defaults: 10,000 independent replications, paired across methods, horizons 50,000, 200,000, 500,000 and 1,000,000, seed 20261018. The polynomial exponent 0.6 satisfies the PR remainder condition with SGD exponent 0.65.

The final three-panel histogram uses common bins and axes and overlays the theoretical normal densities. Raw draws and moment summaries are saved alongside `attainment_histograms.pdf` and `.png`. This experiment illustrates that the estimator with slower single-step MSE can attain the minimum downstream asymptotic variance within the normalized historical linear class.