2018-10-03 18:59:01 +00:00
|
|
|
# setup: import numpy as np ; N = 500000 ; X, Y = np.random.rand(N), np.random.rand(N) # noqa
|
|
|
|
|
# run: lstsqr(X, Y)
|
|
|
|
|
# from:
|
|
|
|
|
# http://nbviewer.ipython.org/github/rasbt/One-Python-benchmark-per-day/blob/master/ipython_nbs/day10_fortran_lstsqr.ipynb
|
2018-04-05 22:07:33 +00:00
|
|
|
|
2018-10-03 18:59:01 +00:00
|
|
|
# pythran export lstsqr(float[], float[])
|
2018-04-05 22:07:33 +00:00
|
|
|
import numpy as np
|
2018-10-03 12:38:48 +00:00
|
|
|
|
|
|
|
|
|
2018-04-05 22:07:33 +00:00
|
|
|
def lstsqr(x, y):
|
|
|
|
|
""" Computes the least-squares solution to a linear matrix equation. """
|
|
|
|
|
x_avg = np.average(x)
|
|
|
|
|
y_avg = np.average(y)
|
|
|
|
|
dx = x - x_avg
|
2020-06-28 18:24:40 +00:00
|
|
|
var_x = np.sum(dx ** 2)
|
2018-04-05 22:07:33 +00:00
|
|
|
cov_xy = np.sum(dx * (y - y_avg))
|
|
|
|
|
slope = cov_xy / var_x
|
2018-10-03 18:59:01 +00:00
|
|
|
y_interc = y_avg - slope * x_avg
|
2018-04-05 22:07:33 +00:00
|
|
|
return (slope, y_interc)
|
