polykin.math

confidence_ellipse

confidence_ellipse(
    ax: Axes,
    center: tuple[float, float],
    cov: Float2x2Matrix,
    ndata: int,
    alpha: float = 0.05,
    color: str = "black",
    label: str | None = None,
    confint: bool = True,
) -> None

Generate a confidence ellipse for 2 jointly estimated parameters using a linear approximation method.

The joint \(100(1-\alpha)\%\) confidence ellipse for the parameters \(\beta=(\beta_1, \beta_2)\) is represented by the domain of values that satisfy the following condition:

\[ \left \{\beta: (\beta -\hat{\beta})^T \hat{V}^{-1}(\beta -\hat{\beta}) \leq 2 F_{2,n-2,1-\alpha} \right \} \]

where \(\hat{\beta}\) is the point estimate of \(\beta\) (obtained by least-squares fitting), \(\hat{V}\) is the scaled variance-covariance matrix, \(n>2\) is the number of data points considered in the regression, \(\alpha\) is the significance level, and \(F\) is the Fisher-Snedecor distribution.

Note

This method is only exact for models that are linear in the parameters. For models that are non-linear in the parameters, the size and shape of the ellipse is only an approximation to the true joint confidence region.

Additionally, the confidence intervals for the individual coefficients are given by:

\[ \hat{\beta}_i \pm \hat{V}_{ii}^{1/2} t_{n-2,1-\alpha/2} \]

where \(t\) is Student's \(t\) distribution.

Reference

• Vugrin, K. W., L. P. Swiler, R. M. Roberts, N. J. Stucky-Mack, and S. P. Sullivan (2007), Confidence region estimation techniques for nonlinear regression in groundwater flow: Three case studies, Water Resour. Res., 43.

PARAMETER	DESCRIPTION
ax	Matplotlib Axes object to which ellipse will be patched. TYPE: Axes
center	Point estimate of the model parameters, \(\hat{\beta}\). TYPE: tuple[float, float]
cov	Scaled variance-covariance matrix, \(\hat{V}\). TYPE: Float2x2Matrix
ndata	Number of data points. TYPE: int
alpha	Significance level, \(\alpha\). TYPE: float DEFAULT: 0.05
color	Color of ellipse contour and center. TYPE: str DEFAULT: 'black'
label	Ellipse label. TYPE: str | None DEFAULT: None
confint	If True the individual confidence intervals of the parameters are represented as lines. TYPE: bool DEFAULT: True

See also

• confidence_region: alternative method based on a rigorous approach for non-linear models.

Examples:

>>> from polykin.math.jcr import confidence_ellipse
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> fig, ax = plt.subplots()
>>> confidence_ellipse(ax=ax, center=(0.3,0.8),
...     cov=np.array([[2e-4, 3e-4], [3e-4, 2e-3]]), ndata=9)

Source code in src/polykin/math/jcr.py

def confidence_ellipse(ax: Axes,
                       center: tuple[float, float],
                       cov: Float2x2Matrix,
                       ndata: int,
                       alpha: float = 0.05,
                       color: str = 'black',
                       label: str | None = None,
                       confint: bool = True
                       ) -> None:
    r"""Generate a confidence ellipse for 2 jointly estimated parameters
    using a linear approximation method.

    The joint $100(1-\alpha)\%$ confidence ellipse for the parameters
    $\beta=(\beta_1, \beta_2)$ is represented by the domain of values that
    satisfy the following condition:

    $$ \left \{\beta: (\beta -\hat{\beta})^T \hat{V}^{-1}(\beta -\hat{\beta})
       \leq 2 F_{2,n-2,1-\alpha} \right \} $$

    where $\hat{\beta}$ is the point estimate of $\beta$ (obtained by
    least-squares fitting), $\hat{V}$ is the _scaled_ variance-covariance
    matrix, $n>2$ is the number of data points considered in the regression,
    $\alpha$ is the significance level, and $F$ is the Fisher-Snedecor
    distribution.

    !!! note

        This method is only exact for models that are linear in the parameters.
        For models that are non-linear in the parameters, the size and shape of
        the ellipse is only an approximation to the true joint confidence
        region.

    Additionally, the confidence intervals for the individual coefficients are
    given by:

    $$ \hat{\beta}_i \pm \hat{V}_{ii}^{1/2} t_{n-2,1-\alpha/2} $$

    where $t$ is Student's $t$ distribution.

    **Reference**

    *   Vugrin, K. W., L. P. Swiler, R. M. Roberts, N. J. Stucky-Mack, and
    S. P. Sullivan (2007), Confidence region estimation techniques for
    nonlinear regression in groundwater flow: Three case studies, Water
    Resour. Res., 43.

    Parameters
    ----------
    ax : Axes
        Matplotlib Axes object to which ellipse will be patched.
    center : tuple[float, float]
        Point estimate of the model parameters, $\hat{\beta}$.
    cov : Float2x2Matrix
        Scaled variance-covariance matrix, $\hat{V}$.
    ndata : int
        Number of data points.
    alpha : float
        Significance level, $\alpha$.
    color : str
        Color of ellipse contour and center.
    label : str | None
        Ellipse label.
    confint : bool
        If `True` the _individual_ confidence intervals of the parameters are
        represented as lines.

    See also
    --------
    * [`confidence_region`](confidence_region.md): alternative method based on
    a rigorous approach for non-linear models.

    Examples
    --------
    >>> from polykin.math.jcr import confidence_ellipse
    >>> import matplotlib.pyplot as plt
    >>> import numpy as np
    >>> fig, ax = plt.subplots()
    >>> confidence_ellipse(ax=ax, center=(0.3,0.8),
    ...     cov=np.array([[2e-4, 3e-4], [3e-4, 2e-3]]), ndata=9)
    """

    # Method implementation is specific for 2D
    npar = 2
    if len(center) != npar:
        raise ShapeError(f"`center` must be a vector of length {npar}.")

    if cov.shape != (npar, npar):
        raise ShapeError(f"`cov` must be a {npar}x{npar} matrix.")

    # Check inputs
    check_bounds(alpha, 0.001, 0.90, 'alpha')
    check_bounds(ndata, npar + 1, np.inf, 'ndata')

    # Eigenvalues and (orthogonal) eigenvectors of cov
    eigen = np.linalg.eigh(cov)
    eigenvector = eigen.eigenvectors[0]
    eigenvalues = eigen.eigenvalues

    # Angle of ellipse wrt to x-axis
    angle = np.degrees(arctan(eigenvector[1]/eigenvector[0]))

    # "Scale" of ellipse
    scale = npar*Fdist.ppf(1. - alpha, npar, ndata - npar)

    # Lengths of ellipse axes
    width, height = 2*sqrt(scale*eigenvalues)

    # Ellipse
    ellipse = Ellipse(center,
                      width=width,
                      height=height,
                      angle=angle,
                      facecolor='none',
                      edgecolor=color,
                      label=label)

    ax.add_patch(ellipse)
    ax.scatter(*center, c=color, s=5)

    # Individual confidence intervals
    if confint:
        ci = sqrt(np.diag(cov))*tdist.ppf(1. - alpha/2, ndata - npar)
        ax.errorbar(*center, xerr=ci[0], yerr=ci[1], color=color)

    return None