0

I'm trying to perform generalized least-squares fitting of highly correlated data. The covariance matrix for these data is large, with many elements close to 1.

# X and Y correspond to my predictor and dependent variables
# C is the covariance matrix of Y
# BF are the best-fit parameters
# BFCM is the covariance matrix of BF

from pylab import *

BF = inv(X.T @ inv(C) @ X) @ (X.T @ inv(C) @ Y)
BFCM = inv(X.T @ inv(C) @ X)

I'm getting negative diagonal terms in the correlation matrix of the best-fit parameters (BFCM), which is definitely wrong. It seems likely that these negative values come from the limited float precision used by numpy.linalg. Is there a way to avoid such floating point precision problems in this case?

Thanks for any suggestions.

Mathieu
  • 241
  • 1
  • 8
  • 2
    Use `numpy.linalg.solve` to start with instead of inverses. – percusse Mar 17 '18 at 11:34
  • Thanks, I'll try that as soon as I can. May I ask why you're suggesting it? – Mathieu Mar 17 '18 at 11:42
  • @Mathieu, https://stackoverflow.com/questions/31256252/why-does-numpy-linalg-solve-offer-more-precise-matrix-inversions-than-numpy-li – MaxU - stand with Ukraine Mar 17 '18 at 11:51
  • @percusse, I tried using `invC = solve(C, eye(C.shape[0]))`; `BFCM = solve(X.T @ invC @ X, eye(2))`; `BF = BFCM @ X.T @ invC @ Y`, but I still get rubbish values for the best-fit parameters. Is this what you were suggesting, or did I misunderstand? – Mathieu Mar 17 '18 at 12:55

0 Answers0