Source code for astroML.linear_model.TLS

import numpy as np


def TLS_logL(v, X, dX):
    """Compute the total least squares log-likelihood

    This uses Hogg et al eq. 29-32

    Parameters
    ----------
    v : ndarray
        The normal vector to the linear best fit.  shape=(D,).
        Note that the magnitude |v| is a stand-in for the intercept.
    X : ndarray
        The input data.  shape = [N, D]
    dX : ndarray
        The covariance of the errors for each point.
        For diagonal errors, the shape = (N, D) and the entries are
        dX[i] = [sigma_x1, sigma_x2 ... sigma_xD]
        For full covariance, the shape = (N, D, D) and the entries are
        dX[i] = Cov(X[i], X[i]), the full error covariance.

    Returns
    -------
    logL : float
        The log-likelihood of the model v given the data.

    Notes
    -----
    This implementation follows Hogg 2010, arXiv 1008.4686
    """
    # check inputs
    X, dX, v = map(np.asarray, (X, dX, v))
    N, D = X.shape
    assert v.shape == (D,)
    assert dX.shape in ((N, D), (N, D, D))

    v_norm = np.linalg.norm(v)
    v_hat = v / v_norm

    # eq. 30
    Delta = np.dot(X, v_hat) - v_norm

    # eq. 31
    if dX.ndim == 2:
        # diagonal covariance
        Sig2 = np.sum(dX * v_hat ** 2, 1)
    else:
        # full covariance
        Sig2 = np.dot(np.dot(v_hat, dX), v_hat)

    return (-0.5 * np.sum(np.log(2 * np.pi * Sig2))
            - np.sum(0.5 * Delta ** 2 / Sig2))