sp_fit_report.py in Python

Fit a learning curve, bootstrap its parameters, and test the residuals.

import numpy as np
from scipy import optimize, stats

CEILING, RATE, NOISE = 118.0, 0.11, 3.5


def learning(x, ceiling, rate):
    """Tokens per minute after x sessions: exponential approach."""
    return ceiling * (1.0 - np.exp(-rate * x))


def make_observations(sessions, rng):
    return learning(sessions, CEILING, RATE) + rng.normal(0, NOISE,
                                                          size=sessions.size)


def fit(sessions, observed):
    params, covariance = optimize.curve_fit(
        learning, sessions, observed, p0=[100.0, 0.05])
    return params, np.sqrt(np.diag(covariance))


def bootstrap_params(sessions, observed, rng, draws=300):
    keep = np.empty((draws, 2))
    n = sessions.size
    for i in range(draws):
        idx = rng.integers(0, n, size=n)
        try:
            keep[i] = fit(sessions[idx], observed[idx])[0]
        except RuntimeError:
            keep[i] = np.nan
    return keep[~np.isnan(keep).any(axis=1)]


def main():
    rng = np.random.default_rng(21)
    sessions = np.arange(1, 61, dtype=float)
    observed = make_observations(sessions, rng)

    params, errors = fit(sessions, observed)
    print("curve_fit")
    print(f"  ceiling {params[0]:7.3f} +/- {errors[0]:.3f}  (true {CEILING})")
    print(f"  rate    {params[1]:7.4f} +/- {errors[1]:.4f}  (true {RATE})")

    draws = bootstrap_params(sessions, observed, rng)
    low, high = np.percentile(draws[:, 0], [2.5, 97.5])
    print(f"bootstrap ceiling 95% interval: {low:.2f} to {high:.2f}"
          f"  ({draws.shape[0]} usable draws)")

    residuals = observed - learning(sessions, *params)
    print("residuals")
    print(f"  sd {residuals.std(ddof=2):.3f} (noise was {NOISE})")
    print(f"  shapiro p {stats.shapiro(residuals).pvalue:.4f}")
    print(f"  runs correlation {stats.pearsonr(sessions, residuals)[1]:.4f}")

    at_90 = optimize.root_scalar(
        lambda x: learning(x, *params) - 0.9 * params[0],
        bracket=[1, 500], method="brentq")
    print(f"90% of the ceiling reached by session {at_90.root:.1f}")

    total = float(np.trapezoid(learning(sessions, *params), sessions))
    print(f"area under the fitted curve: {total:.0f}")


if __name__ == "__main__":
    main()

How it works

1. curve_fit gives the parameters; the covariance gives their error.
2. A bootstrap over resampled rows checks that error empirically.
3. The residuals are then tested for normality and structure.

Keywords and builtins used here

• as
• bootstrap_params
• def
• except
• fit
• float
• for
• if
• lambda
• learning
• main
• make_observations
• print
• range
• return
• try

The run, in numbers

Lines

66

Characters to type

1964

Tokens

576

Three-star pace

115 tpm

At the three-star pace of 115 tokens a minute, this run takes about 301 seconds.

Type this snippet

Step 1 of 1 in Encore, step 23 of 23 in Scientific computing with SciPy.

← Previous