Skip to content

Bootstrap

semiparametric_bootstrap()

Semiparametric bootstrap for parameter uncertainty.

Resamples death counts from their assumed distribution (Poisson or Binomial) conditional on the fitted expected deaths, then refits the model on each bootstrap sample.

Parameters:

Name Type Description Default
fit FitStMoMo

Fitted model to bootstrap.

 required
nboot int

Number of bootstrap replicates.

 500
seed int | None

Random seed.

 None
n_jobs int

Number of parallel jobs. Requires joblib (pip install pystmomo[parallel]). Defaults to 1 (sequential).

 1

Returns:

Type Description
BootStMoMo

Examples:

>>> from pystmomo import lc, load_ew_male, semiparametric_bootstrap
>>> data = load_ew_male()
>>> fit = lc().fit(data.deaths, data.exposures, ages=data.ages, years=data.years)
>>> boot = semiparametric_bootstrap(fit, nboot=100, seed=0)
>>> lo, hi = boot.parameter_ci("kt", level=0.95)
>>> lo.shape
(1, 51)
Source code in src/pystmomo/bootstrap/semipar_boot.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def semiparametric_bootstrap(
    fit: FitStMoMo,
    nboot: int = 500,
    *,
    seed: int | None = None,
    n_jobs: int = 1,
) -> BootStMoMo:
    """Semiparametric bootstrap for parameter uncertainty.

    Resamples death counts from their assumed distribution (Poisson or
    Binomial) conditional on the fitted expected deaths, then refits the
    model on each bootstrap sample.

    Parameters
    ----------
    fit:
        Fitted model to bootstrap.
    nboot:
        Number of bootstrap replicates.
    seed:
        Random seed.
    n_jobs:
        Number of parallel jobs.  Requires ``joblib`` (``pip install
        pystmomo[parallel]``).  Defaults to 1 (sequential).

    Returns
    -------
    BootStMoMo

    Examples
    --------
    >>> from pystmomo import lc, load_ew_male, semiparametric_bootstrap
    >>> data = load_ew_male()
    >>> fit = lc().fit(data.deaths, data.exposures, ages=data.ages, years=data.years)
    >>> boot = semiparametric_bootstrap(fit, nboot=100, seed=0)
    >>> lo, hi = boot.parameter_ci("kt", level=0.95)
    >>> lo.shape
    (1, 51)
    """
    rng = np.random.default_rng(seed)
    seeds = rng.integers(0, 2**31, size=nboot).tolist()

    def _one_replicate(s: int) -> FitStMoMo | None:
        rng_b = np.random.default_rng(s)
        try:
            D_boot = _sample_deaths(fit, rng_b)
            return fit.model.fit(
                D_boot, fit.Ext, fit.ages, fit.years,
                wxt=fit.wxt, oxt=fit.oxt,
            )
        except Exception:
            return None

    if n_jobs != 1:
        try:
            from joblib import Parallel, delayed
            results = Parallel(n_jobs=n_jobs)(
                delayed(_one_replicate)(s) for s in seeds
            )
        except ImportError:
            results = [_one_replicate(s) for s in seeds]
    else:
        results = [_one_replicate(s) for s in seeds]

    valid_fits = [r for r in results if r is not None]
    boot = BootStMoMo(
        base_fit=fit,
        nboot=nboot,
        method="semiparametric",
        fits=valid_fits,
    )
    return boot

residual_bootstrap()

Residual bootstrap for parameter uncertainty.

Resamples deviance residuals from the fitted model, converts them back to death counts, and refits the model on each bootstrap sample.

Parameters:

Name Type Description Default
fit FitStMoMo

Fitted model to bootstrap.

 required
nboot int

Number of bootstrap replicates.

 500
seed int | None

Random seed.

 None
n_jobs int

Number of parallel jobs (requires joblib).

 1

Returns:

Type Description
BootStMoMo
Source code in src/pystmomo/bootstrap/residual_boot.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def residual_bootstrap(
    fit: FitStMoMo,
    nboot: int = 500,
    *,
    seed: int | None = None,
    n_jobs: int = 1,
) -> BootStMoMo:
    """Residual bootstrap for parameter uncertainty.

    Resamples deviance residuals from the fitted model, converts them back to
    death counts, and refits the model on each bootstrap sample.

    Parameters
    ----------
    fit:
        Fitted model to bootstrap.
    nboot:
        Number of bootstrap replicates.
    seed:
        Random seed.
    n_jobs:
        Number of parallel jobs (requires ``joblib``).

    Returns
    -------
    BootStMoMo
    """
    rng = np.random.default_rng(seed)
    seeds = rng.integers(0, 2**31, size=nboot).tolist()

    # Compute deviance residuals once
    res = deviance_residuals(fit)   # (n_ages, n_years), masked cells = 0

    def _one_replicate(s: int) -> FitStMoMo | None:
        rng_b = np.random.default_rng(s)
        try:
            D_boot = _resample_residuals(fit, res, rng_b)
            return fit.model.fit(
                D_boot, fit.Ext, fit.ages, fit.years,
                wxt=fit.wxt, oxt=fit.oxt,
            )
        except Exception:
            return None

    if n_jobs != 1:
        try:
            from joblib import Parallel, delayed
            results = Parallel(n_jobs=n_jobs)(
                delayed(_one_replicate)(s) for s in seeds
            )
        except ImportError:
            results = [_one_replicate(s) for s in seeds]
    else:
        results = [_one_replicate(s) for s in seeds]

    valid_fits = [r for r in results if r is not None]
    return BootStMoMo(
        base_fit=fit,
        nboot=nboot,
        method="residual",
        fits=valid_fits,
    )

BootStMoMo

Result of bootstrap uncertainty quantification.

Attributes:

Name Type Description
base_fit FitStMoMo

The original fitted model.
nboot int

Number of bootstrap replicates.
method Literal['semiparametric', 'residual']

Bootstrap method used: "semiparametric" or "residual".
fits list[FitStMoMo]

List of refitted :class:FitStMoMo objects (one per replicate).
ax_b ndarray | None

Bootstrap distribution of α_x, shape (nboot, n_ages) or None.
bx_b ndarray

Bootstrap distribution of β_x, shape (nboot, n_ages, N).
kt_b ndarray

Bootstrap distribution of κ_t, shape (nboot, N, n_years).
b0x_b ndarray | None

Bootstrap distribution of β_x^(0), shape (nboot, n_ages) or None.
gc_b ndarray | None

Bootstrap distribution of γ_c, shape (nboot, n_cohorts) or None.
Source code in src/pystmomo/bootstrap/boot_result.py 
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
@dataclass
class BootStMoMo:
    """Result of bootstrap uncertainty quantification.

    Attributes
    ----------
    base_fit:
        The original fitted model.
    nboot:
        Number of bootstrap replicates.
    method:
        Bootstrap method used: ``"semiparametric"`` or ``"residual"``.
    fits:
        List of refitted :class:`FitStMoMo` objects (one per replicate).
    ax_b:
        Bootstrap distribution of α_x, shape (nboot, n_ages) or None.
    bx_b:
        Bootstrap distribution of β_x, shape (nboot, n_ages, N).
    kt_b:
        Bootstrap distribution of κ_t, shape (nboot, N, n_years).
    b0x_b:
        Bootstrap distribution of β_x^(0), shape (nboot, n_ages) or None.
    gc_b:
        Bootstrap distribution of γ_c, shape (nboot, n_cohorts) or None.
    """

    base_fit: FitStMoMo
    nboot: int
    method: Literal["semiparametric", "residual"]
    fits: list[FitStMoMo] = field(default_factory=list)

    @property
    def ax_b(self) -> np.ndarray | None:
        """Bootstrap ax, shape (nboot, n_ages)."""
        if not self.fits or self.fits[0].ax is None:
            return None
        return np.stack([f.ax for f in self.fits])

    @property
    def bx_b(self) -> np.ndarray:
        """Bootstrap bx, shape (nboot, n_ages, N)."""
        return np.stack([f.bx for f in self.fits])

    @property
    def kt_b(self) -> np.ndarray:
        """Bootstrap kt, shape (nboot, N, n_years)."""
        return np.stack([f.kt for f in self.fits])

    @property
    def b0x_b(self) -> np.ndarray | None:
        """Bootstrap b0x, shape (nboot, n_ages) or None."""
        if not self.fits or self.fits[0].b0x is None:
            return None
        return np.stack([f.b0x for f in self.fits])

    @property
    def gc_b(self) -> np.ndarray | None:
        """Bootstrap gc, shape (nboot, n_cohorts) or None."""
        if not self.fits or self.fits[0].gc is None:
            return None
        return np.stack([f.gc for f in self.fits])

    def parameter_ci(
        self,
        param: str = "kt",
        level: float = 0.95,
    ) -> tuple[np.ndarray, np.ndarray]:
        """Bootstrap confidence interval for a parameter.

        Parameters
        ----------
        param:
            One of ``"ax"``, ``"bx"``, ``"kt"``, ``"b0x"``, ``"gc"``.
        level:
            Confidence level in (0, 1).

        Returns
        -------
        lower, upper:
            Arrays of the same shape as the parameter, at the requested
            confidence level.
        """
        alpha = (1.0 - level) / 2.0
        arr = getattr(self, f"{param}_b")
        if arr is None:
            raise ValueError(f"Parameter '{param}' is not available in this bootstrap.")
        lower = np.quantile(arr, alpha, axis=0)
        upper = np.quantile(arr, 1 - alpha, axis=0)
        return lower, upper

    def __repr__(self) -> str:
        return (
            f"BootStMoMo(method={self.method!r}, nboot={self.nboot}, "
            f"n_fits_ok={len(self.fits)})"
        )

ax_b property

Bootstrap ax, shape (nboot, n_ages).

b0x_b property

Bootstrap b0x, shape (nboot, n_ages) or None.

bx_b property

Bootstrap bx, shape (nboot, n_ages, N).

gc_b property

Bootstrap gc, shape (nboot, n_cohorts) or None.

kt_b property

Bootstrap kt, shape (nboot, N, n_years).

parameter_ci(param='kt', level=0.95)

Bootstrap confidence interval for a parameter.

Parameters:

Name Type Description Default
param str

One of "ax", "bx", "kt", "b0x", "gc".

 'kt'
level float

Confidence level in (0, 1).

 0.95

Returns:

Type Description
lower, upper:

Arrays of the same shape as the parameter, at the requested confidence level.
Source code in src/pystmomo/bootstrap/boot_result.py 
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def parameter_ci(
    self,
    param: str = "kt",
    level: float = 0.95,
) -> tuple[np.ndarray, np.ndarray]:
    """Bootstrap confidence interval for a parameter.

    Parameters
    ----------
    param:
        One of ``"ax"``, ``"bx"``, ``"kt"``, ``"b0x"``, ``"gc"``.
    level:
        Confidence level in (0, 1).

    Returns
    -------
    lower, upper:
        Arrays of the same shape as the parameter, at the requested
        confidence level.
    """
    alpha = (1.0 - level) / 2.0
    arr = getattr(self, f"{param}_b")
    if arr is None:
        raise ValueError(f"Parameter '{param}' is not available in this bootstrap.")
    lower = np.quantile(arr, alpha, axis=0)
    upper = np.quantile(arr, 1 - alpha, axis=0)
    return lower, upper