Skip to content

gelman_rubin #

Compute the Gelman-Rubin diagnostic.

Functions:

gelman_rubin #

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

Compute the Gelman-Rubin diagnostic according to equation 4 in Statist. Sci. 36(4): 518-529 (November 2021). DOI: 10.1214/20-STS812

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 Gelman-Rubin diagnostic.

Source code in red/gelman_rubin.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
def gelman_rubin(data: _npt.NDArray[_np.float64]) -> float:
    """
    Compute the Gelman-Rubin diagnostic according to
    equation 4 in  Statist. Sci. 36(4): 518-529
    (November 2021). DOI: 10.1214/20-STS812

    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 Gelman-Rubin diagnostic.
    """
    # Check that the data is valid.
    data = _check_data(data, one_dim_allowed=False)
    _, n_samples = data.shape

    # Compute the variance estimates.
    intra_run_variance_est = _intra_run_variance(data)
    inter_run_variance_est = _inter_run_variance(data)

    # Get combined variance estimate.
    combined_variance_est = (
        n_samples - 1
    ) / n_samples * intra_run_variance_est + inter_run_variance_est / n_samples

    # Compute GR diagnostic.
    gelman_rubin_diagnostic = _np.sqrt(combined_variance_est / intra_run_variance_est)

    return gelman_rubin_diagnostic  # type: ignore[no-any-return]

stable_gelman_rubin #

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

Compute the stable Gelman-Rubin diagnostic according to equation 7 in Statist. Sci. 36(4): 518-529 (November 2021). DOI: 10.1214/20-STS812. This is applicable to a single run.

Source code in red/gelman_rubin.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def stable_gelman_rubin(data: _npt.NDArray[_np.float64], n_pow: float = 1 / 3) -> float:
    """
    Compute the stable Gelman-Rubin diagnostic according to
    equation 7 in  Statist. Sci. 36(4): 518-529
    (November 2021). DOI: 10.1214/20-STS812. This is applicable to
    a single run.
    """
    # Validate the data.
    data = _check_data(data, one_dim_allowed=True)
    _, n_samples = data.shape

    # Compute the variance estimates.
    intra_run_variance_est = _intra_run_variance(data)
    lugsail_variance_est = _lugsail_variance(data, n_pow=n_pow)

    # Get combined variance estimate.
    combined_variance_est = (
        n_samples - 1
    ) / n_samples * intra_run_variance_est + lugsail_variance_est / n_samples

    # Compute GR diagnostic.
    gelman_rubin_diagnostic = _np.sqrt(combined_variance_est / intra_run_variance_est)

    return gelman_rubin_diagnostic  # type: ignore[no-any-return]