Viděl jsem tento seznam zde a nemohl jsem uvěřit, že existuje tolik způsobů řešení nejmenších čtverců. „Normální rovnice“ na Wikipedii se zdály být docela přímočarým způsobem: $$ {\ displaystyle {\ begin {aligned} {\ hat {\ alpha}} & = {\ bar {y}} - {\ hat {\ beta}} \, {\ bar {x}}, \\ {\ hat {\ beta}} & = {\ frac {\ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}}) (y_ {i} - {\ bar {y}} )} {\ sum _ {i = 1} ^ {n} (x_ {i} - {\ bar {x}}) ^ {2}}} \ end {zarovnáno}}} $$
Proč je tedy nepoužívat? Předpokládal jsem, že musí existovat problém s výpočtem nebo přesností, protože v prvním odkazu výše Mark L. Stone uvádí, že SVD nebo QR jsou populární metody statistického softwaru a že normální rovnice jsou „TERRIBLE z hlediska spolehlivosti a numerické přesnosti“. Nicméně v následujícím kódu mi normální rovnice dávají přesnost na ~ 12 desetinných míst ve srovnání se třemi populárními funkcemi pythonu: numpy's polyfit; linregress společnosti scipy; a LinearRegression společnosti scikit-learn.
Zajímavější je, že metoda běžných rovnic je nejrychlejší, když n = 100000000. Výpočetní časy pro mě jsou: 2,5 s pro linregress; 12,9 s pro polyfit; 4,2 s pro lineární regrese; a 1,8 s pro normální rovnici.
Kód:
import numpy jako np
ze sklearn.linear_model import LinearRegression
ze scipy.stats import linregress
importovat čas
b0 = 0
b1 = 1
n = 100000000
x = np.linspace (-5, 5, n)
np.random.seed (42)
e = np.random.randn (n)
y = b0 + b1 * x + e
# scipy
start = timeit.default_timer ()
print (str.format ('{0: .30f}', linregress (x, y) [0])))
stop = timeit.default_timer ()
tisk (stop - start)
# numpy
start = timeit.default_timer ()
print (str.format ('{0: .30f}', np.polyfit (x, y, 1) [0]))
stop = timeit.default_timer ()
tisk (stop - start)
# sklearn
clf = LinearRegression ()
start = timeit.default_timer ()
clf.fit (x.reshape (-1, 1), y.reshape (-1, 1))
stop = timeit.default_timer ()
print (str.format ('{0: .30f}', clf.coef_ [0, 0]))
tisk (stop - start)
# normální rovnice
start = timeit.default_timer ()
sklon = np.sum ((x-x.mean ()) * (y-y.mean ())) / np.sum ((x-x.mean ()) ** 2)
stop = timeit.default_timer ()
print (str.format ('{0: .30f}', sklon))
tisk (stop - start)