Skip to content

ess #

Functions to calculate the statistical inefficiency and effective sample size.

Functions:

convert_sse_series_to_ess_series #

convert_sse_series_to_ess_series(
    data: NDArray[float64], sse_series: NDArray[float64]
) -> NDArray[float64]

Convert a series of squared standard errors to a series of effective sample sizes.

Parameters:

  • data (ndarray) –

    The time series data with shape (n_chains, n_samples).

  • sse_series (ndarray) –

    The squared standard error series.

Returns:

  • ndarray –

    The effective sample size series.

Source code in red/ess.py 
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
def convert_sse_series_to_ess_series(
    data: _npt.NDArray[_np.float64], sse_series: _npt.NDArray[_np.float64]
) -> _npt.NDArray[_np.float64]:
    """Convert a series of squared standard errors to a series of effective sample sizes.

    Parameters
    ----------
    data : np.ndarray
        The time series data with shape (n_chains, n_samples).

    sse_series : np.ndarray
        The squared standard error series.

    Returns
    -------
    np.ndarray
        The effective sample size series.

    """
    # Validate the data.
    data = check_data(data, one_dim_allowed=True)

    # Now get the uncorrelated variances.
    uncor_vars = _np.zeros_like(sse_series)

    for i in range(len(sse_series)):
        # Get "biased", rather than n - 1, variance.
        uncor_vars[i] = data[:, i:].var()

    return uncor_vars / sse_series

get_ess_series_init_seq #

get_ess_series_init_seq(
    data: NDArray[float64],
    sequence_estimator: str = "initial_convex",
    min_max_lag_time: int = 3,
    max_max_lag_time: int | None = None,
    smooth_lag_times: bool = False,
    frac_padding: float = 0.1,
) -> tuple[NDArray[float64], NDArray[float64]]

Compute a series of effective sample sizes for a time series.

Effective sample sizes are computed as data is discarded from the beginning of the time series. The autocorrelation is computed using the sequence estimator specified.

Parameters:

  • data (ndarray) –

    A time series of data with shape (n_samples,).

  • sequence_estimator (str, default: 'initial_convex' ) –

    The initial sequence estimator to use. Can be "positive", "initial_positive", "initial_monotone", or "initial_convex". The default is "initial_convex". "positive" corresponds to truncating the auto-covariance function at the first negative value, as is done in pymbar. The other methods correspond to the methods described in Geyer, 1992: https://www.jstor.org/stable/2246094.

  • min_max_lag_time (int, default: 3 ) –

    The minimum maximum lag time to use when estimating the statistical inefficiency.

  • max_max_lag_time (int, default: None ) –

    The maximum maximum lag time to use when calculating the auto-correlation function. If None, the maximum lag time will be the length of the time series.

  • smooth_lag_times (bool, default: False ) –

    Whether to smooth out the max lag times by a) converting them to a monotinically decreasing sequence and b) linearly interpolating between points where the sequence changes. This may be useful when the max lag times are noisy.

  • frac_padding (float, default: 0.1 ) –

    The fraction of the end of the timeseries to avoid calculating the variance for. For example, if frac_padding = 0.1, the variance will be calculated for the first 90% of the time series. This helps to avoid noise in the variance when there are few data points.

Returns:

  • ndarray –

    The effective sample size series.

  • ndarray –

    The maximum lag times used.

Source code in red/ess.py 
 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
108
109
def get_ess_series_init_seq(
    data: _npt.NDArray[_np.float64],
    sequence_estimator: str = "initial_convex",
    min_max_lag_time: int = 3,
    max_max_lag_time: int | None = None,
    smooth_lag_times: bool = False,
    frac_padding: float = 0.1,
) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
    """Compute a series of effective sample sizes for a time series.

    Effective sample sizes are computed as data is discarded from the beginning of the
    time series. The autocorrelation is computed using the sequence estimator specified.

    Parameters
    ----------
    data : numpy.ndarray
        A time series of data with shape (n_samples,).

    sequence_estimator : str, optional
        The initial sequence estimator to use. Can be "positive", "initial_positive",
        "initial_monotone", or "initial_convex". The default is "initial_convex". "positive"
        corresponds to truncating the auto-covariance function at the first negative value, as is
        done in pymbar. The other methods correspond to the methods described in Geyer, 1992:
        https://www.jstor.org/stable/2246094.

    min_max_lag_time : int, optional, default=3
        The minimum maximum lag time to use when estimating the statistical inefficiency.

    max_max_lag_time : int, optional, default=None
        The maximum maximum lag time to use when calculating the auto-correlation function.
        If None, the maximum lag time will be the length of the time series.

    smooth_lag_times : bool, optional, default=False
        Whether to smooth out the max lag times by a) converting them to a monotinically
        decreasing sequence and b) linearly interpolating between points where the sequence
        changes. This may be useful when the max lag times are noisy.

    frac_padding : float, optional, default=0.1
        The fraction of the end of the timeseries to avoid calculating the variance
        for. For example, if frac_padding = 0.1, the variance will be calculated
        for the first 90% of the time series. This helps to avoid noise in the
        variance when there are few data points.

    Returns
    -------
    np.ndarray
        The effective sample size series.

    np.ndarray
        The maximum lag times used.

    """
    sse_series, max_lag_times = _get_sse_series_init_seq(
        data,
        sequence_estimator=sequence_estimator,
        min_max_lag_time=min_max_lag_time,
        max_max_lag_time=max_max_lag_time,
        smooth_lag_times=smooth_lag_times,
        frac_padding=frac_padding,
    )

    ess_series = convert_sse_series_to_ess_series(data, sse_series)

    return ess_series, max_lag_times

get_ess_series_window #

get_ess_series_window(
    data: NDArray[float64],
    kernel: Callable[[int], NDArray[float64]] = bartlett,
    window_size_fn: (
        Callable[[int], int] | None
    ) = lambda x: round(x**0.5),
    window_size: int | None = None,
) -> tuple[NDArray[float64], NDArray[float64]]

Compute a series of effective sample sizes for a time series.

Effective sample sizes are computed as data is discarded from the beginning of the time series. The squared standard error is computed using the window size and kernel specified.

Parameters:

  • data (ndarray) –

    A time series of data with shape (n_samples,).

  • kernel (callable, default: numpy.bartlett ) –

    A function that takes a window size and returns a window function.

  • window_size_fn (callable, default: lambda x: round(x**0.5) ) –

    A function that takes the length of the time series and returns the window size to use. If this is not None, window_size must be None.

  • window_size (int, default: None ) –

    The size of the window to use, defined in terms of time lags in the forwards direction. If this is not None, window_size_fn must be None.

Returns:

  • ndarray –

    The squared standard error series.

  • ndarray –

    The window sizes used.

Source code in red/ess.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def get_ess_series_window(
    data: _npt.NDArray[_np.float64],
    kernel: _Callable[[int], _npt.NDArray[_np.float64]] = _np.bartlett,
    window_size_fn: _Callable[[int], int] | None = lambda x: round(x**0.5),
    window_size: int | None = None,
) -> tuple[_npt.NDArray[_np.float64], _npt.NDArray[_np.float64]]:
    """Compute a series of effective sample sizes for a time series.

    Effective sample sizes are computed as data is discarded from the beginning of the
    time series. The squared standard error is computed using the window size and kernel
    specified.

    Parameters
    ----------
    data : numpy.ndarray
        A time series of data with shape (n_samples,).

    kernel : callable, optional, default=numpy.bartlett
        A function that takes a window size and returns a window function.

    window_size_fn : callable, optional, default=lambda x: round(x**0.5)
        A function that takes the length of the time series and returns the window size
        to use. If this is not None, window_size must be None.

    window_size : int, optional, default=None
        The size of the window to use, defined in terms of time lags in the
        forwards direction. If this is not None, window_size_fn must be None.

    Returns
    -------
    np.ndarray
        The squared standard error series.

    np.ndarray
        The window sizes used.

    """
    sse_series, max_lag_times = _get_sse_series_window(
        data, kernel=kernel, window_size_fn=window_size_fn, window_size=window_size
    )

    ess_series = convert_sse_series_to_ess_series(data, sse_series)

    return ess_series, max_lag_times

statistical_inefficiency_inter_variance #

statistical_inefficiency_inter_variance(
    data: NDArray[float64],
) -> float

Compute the statistical inefficiency of a time series.

The statistical inefficiency is computed by dividing the inter-run variance estimate by the intra-run variance estimate. More than one run is required.

Parameters:

  • data (ndarray) –

    The time series data. This should have shape (n_runs, n_samples) and must have at least two runs.

Returns:

  • float –

    The statistical inefficiency.

Source code in red/ess.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def statistical_inefficiency_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
    """Compute the statistical inefficiency of a time series.

    The statistical inefficiency is computed by dividing the inter-run variance estimate
    by the intra-run variance estimate. More than one run is required.

    Parameters
    ----------
    data : np.ndarray
        The time series data. This should have shape
        (n_runs, n_samples) and must have at least two
        runs.

    Returns
    -------
    float
        The statistical inefficiency.

    """
    g = inter_run_variance(data) / intra_run_variance(data)
    # Ensure that the statistical inefficiency is at least 1.
    return max(g, 1)

statistical_inefficiency_lugsail_variance #

statistical_inefficiency_lugsail_variance(
    data: NDArray[float64], n_pow: float = 1 / 3
) -> float

Compute the statistical inefficiency of a time series.

The statistical inefficiency is computed by dividing the lugsail replicated batch means variance estimate by the intra-run variance estimate. This is applicable to a single run.

Parameters:

  • data (ndarray) –

    The time series data. This should have shape (n_runs, n_samples) and must have at least two runs.

  • n_pow (float, default: 1 / 3 ) –

    The power to use in the lugsail variance estimate. This should be between 0 and 1. The default is 1/3.

Returns:

  • float –

    The statistical inefficiency.

Source code in red/ess.py 
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
def statistical_inefficiency_lugsail_variance(
    data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3
) -> float:
    """Compute the statistical inefficiency of a time series.

    The statistical inefficiency is computed by dividing the lugsail replicated batch
    means variance estimate by the intra-run variance estimate. This is applicable to a
    single run.

    Parameters
    ----------
    data : np.ndarray
        The time series data. This should have shape
        (n_runs, n_samples) and must have at least two
        runs.

    n_pow : float, optional
        The power to use in the lugsail variance estimate. This
        should be between 0 and 1. The default is 1/3.

    Returns
    -------
    float
        The statistical inefficiency.

    """
    g = lugsail_variance(data, n_pow=n_pow) / intra_run_variance(data)
    # Ensure that the statistical inefficiency is at least 1.
    return max(g, 1)

ess_inter_variance #

ess_inter_variance(data: NDArray[float64]) -> float

Compute the effective sample size of a time series.

The effective sample size is computed by dividing the total number of samples by the statistical inefficiency, where the statistical inefficiency is calculated using the ratio of the inter-run and intra-run variance estimates.

Parameters:

  • data (ndarray) –

    The time series data. This should have shape (n_runs, n_samples) and must have at least two runs.

Returns:

  • float –

    The effective sample size.

Source code in red/ess.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
def ess_inter_variance(data: _npt.NDArray[_np.float64]) -> float:
    """Compute the effective sample size of a time series.

    The effective sample size is computed by dividing the total number of samples by the
    statistical inefficiency, where the statistical inefficiency is calculated using the
    ratio of the inter-run and intra-run variance estimates.

    Parameters
    ----------
    data : np.ndarray
        The time series data. This should have shape
        (n_runs, n_samples) and must have at least two
        runs.

    Returns
    -------
    float
        The effective sample size.

    """
    data = check_data(data, one_dim_allowed=False)
    n_runs: int = data.shape[0]
    n_samples: int = data.shape[1]
    total_samples = n_runs * n_samples
    return total_samples / statistical_inefficiency_inter_variance(data)

ess_lugsail_variance #

ess_lugsail_variance(
    data: NDArray[float64], n_pow: float = 1 / 3
) -> float

Compute the effective sample size of a time series.

The effective sample size is computed by dividing the total number of samples by the statistical inefficiency, where the statistical inefficiency is calculated using the ratio of the lugsail replicated batch means and intra-run variance estimates.

Parameters:

  • data (ndarray) –

    The time series data. This should have shape (n_runs, n_samples) and must have at least two runs.

  • n_pow (float, default: 1 / 3 ) –

    The power to use in the lugsail variance estimate. This should be between 0 and 1. The default is 1/3.

Returns:

  • float –

    The effective sample size.

Source code in red/ess.py 
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
def ess_lugsail_variance(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> float:
    """Compute the effective sample size of a time series.

    The effective sample size is computed by dividing the total number of samples by the
    statistical inefficiency, where the statistical inefficiency is calculated using the
    ratio of the lugsail replicated batch means and intra-run variance estimates.

    Parameters
    ----------
    data : np.ndarray
        The time series data. This should have shape
        (n_runs, n_samples) and must have at least two
        runs.

    n_pow : float, optional
        The power to use in the lugsail variance estimate. This
        should be between 0 and 1. The default is 1/3.

    Returns
    -------
    float
        The effective sample size.

    """
    data = check_data(data, one_dim_allowed=True)
    n_runs: int = data.shape[0]
    n_samples: int = data.shape[1]
    total_samples = n_runs * n_samples
    return total_samples / statistical_inefficiency_lugsail_variance(data, n_pow=n_pow)