clophfit.fitting.robust_se ========================== .. py:module:: clophfit.fitting.robust_se .. autoapi-nested-parse:: Heteroskedasticity-robust (sandwich) standard errors for a fitted well. The usual standard error trusts the weights: it is ``s^2 (J'J)^-1``, which is only right when the residuals really do have the variance the weights assumed. On these plates they do not - the scatter grows with signal, and with plate position - yet weighting sigma to match costs precision in K, because it down-weights the plateau points that fix S0 and S1 (see :mod:`clophfit.fitting.gain_calibration`). The sandwich estimator takes the other route: keep the weights that fit best and repair only the uncertainty. With ``J`` the weighted Jacobian and ``e`` the weighted residuals, .. math:: \widehat{\operatorname{Var}}(\hat\beta) = (J'J)^{-1} J' \Omega J (J'J)^{-1} where ``Omega`` is diagonal in the observed squared residuals rather than in an assumed sigma: ``e_i^2`` (HC0), ``e_i^2 / (1 - h_i)`` (HC2) or ``e_i^2 / (1 - h_i)^2`` (HC3), ``h`` the leverage. HC3 is the default here because the series are short - seven points and three parameters per label - and HC0 is badly biased in small samples. What this buys: an interval for K that reflects the scatter the well actually shows, including the part the noise model does not describe. What it does not buy: a different K. **Read the size of the correction with care on short series.** On simulated seven-step, two-label wells (nine residual degrees of freedom) HC3 inflates K's standard error by 1.25x whether the noise is homoscedastic or strongly signal-dependent, and it lifts 95% coverage from 0.91 to 0.95 - but so does using a Student-t quantile with the classical error, which covers 0.94-0.95 while HC3 with the same quantile over-covers at 0.96-0.97. At this series length the inflation is mostly a small-sample correction and not evidence of misspecification; the sandwich earns its keep where the design is longer or the wells are pooled. Functions --------- .. autoapisummary:: clophfit.fitting.robust_se.hc_covariance clophfit.fitting.robust_se.sandwich_se clophfit.fitting.robust_se.robust_k_se Module Contents --------------- .. py:function:: hc_covariance(jac, resid, *, kind = 'HC3') Sandwich covariance of the parameters, from the residuals as they are. :param jac: ``(n, p)`` Jacobian of the residuals actually minimised (weighted, if the fit was weighted). :type jac: np.ndarray :param resid: Those residuals, on the same scale. :type resid: np.ndarray :param kind: Small-sample correction: ``HC0`` none, ``HC2`` divides each squared residual by ``1 - h``, ``HC3`` by ``(1 - h)^2`` (the default). :type kind: HCKind :returns: ``(p, p)`` covariance, or None when the Jacobian is rank deficient. :rtype: np.ndarray | None .. rubric:: Examples A straight line whose noise grows with x: the robust slope error exceeds the classical one, which believes the constant variance it was given. >>> import numpy as np >>> rng = np.random.default_rng(0) >>> x = np.linspace(0, 1, 200) >>> j = np.column_stack([np.ones_like(x), x]) >>> e = rng.normal(0, 0.1 + 2 * x) >>> cov = hc_covariance(j, e - j @ np.linalg.lstsq(j, e, rcond=None)[0]) >>> classical = np.var(e, ddof=2) * np.linalg.inv(j.T @ j) >>> bool(cov[1, 1] > classical[1, 1]) True .. py:function:: sandwich_se(jac, resid, *, kind = 'HC3') Robust standard error per parameter (the square root of :func:`hc_covariance`). :param jac: Weighted Jacobian. :type jac: np.ndarray :param resid: Weighted residuals. :type resid: np.ndarray :param kind: Small-sample correction. :type kind: HCKind :returns: Standard errors, or None when the Jacobian is rank deficient. :rtype: np.ndarray | None .. py:function:: robust_k_se(table, params, *, is_ph = True, kind = 'HC3') Model-based and robust standard errors of K, per well. Each well is linearised at its fitted parameters (the same weighted Jacobian the residual tests use, residuals projected through ``I - H`` so a posterior mean is treated as its least-squares equivalent). K is the first column of that design, so its standard error is read from the two covariances: the classical one, ``s^2 (J'J)^-1``, which trusts the weights, and the sandwich one, which does not. :param table: Residual table with ``well``, ``label``, ``x``, ``raw_res``, ``sigma``. :type table: pd.DataFrame :param params: Well to its fitted ``K``, ``S0_