1
2
3
4
5
6
7
8
9
10
11
12
13
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
44
45
46
47
48
49
50
51
52
53
|
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))
|
