Linearna regresija metodom najmanjih kvadrata

Linearna regresija predstavlja linearni pristup modeliranju relacije između jedne ili više nezavisnih varijabli sa zavisnom varijablom. Kada govorimo o zavisnosti varijable $\bf{Y}$ od jedne nezavisne $\bf{X}$ tada govorimo o prostoj linearnoj regresiji. Takav slučaj će biti razmatran u tekstu koji slijedi.

Definirajmo linearnu vezu između veličina $\bf{X}$ i $\bf{Y}$ kao:

$$ Y=a_0+a_1X $$

gdje $a_0$ predstavlja presjek prave $Y$ sa x-osom, a $a_1$ nagib prave. Ovu jednačinu ćemo koristiti pri modeliranju relacije između dvije veličine. Potrebno je odrediti vrijednosti $a_0$ i $a_1$, koristeći dostupne podatke (vrijednosti za $\bf{X}$ i $\bf{Y}$ kako bi dobili pravac $Y(X)$ koji najbolje opisuje relaciju između dvije varijable.

Kako bi minimizirali grešku koju ćemo pri tome napraviti (podaci obično odstupaju od linearne zavisnosti) potrebno je izračunati grešku. Nekoliko je načina na koji se greška može izračunati tako što se definiše odgovarajuća funkcija pogreške. Metoda najmanjih kvadrata podrazumijeva da se funkcija: $$ L=\sum_{i=1}^n (y_i - p_i)^2 $$ gdje su $y_i$ date vrijednosti (npr. eksperimentalni podaci), a $p_i$ vrijednosti koje će dati naš model. Minimiziranjem funkcije (L(x,y)), tj. određivanjem parcijalnih izvoda i njihovim izjednačavanjem sa nulom, dobiju se odgovarajuće vrijednosti za koeficijente $a_0$ i $a_1$ koji će dati najmanju grešku.

Nakon nekoliko linija računa stiže se do relacija za koeficijente: $$ a_1 = \frac{\sum_{i=1}^n (x_i - \bar x)(y_i - \bar y)}{\sum_{i=1}^n (x_i - \bar x)^2}$$

$$ a_0 = \bar y - m\bar x$$

gdje $\bar{x}$ i $\bar{y}$ predstavljaju odgovarajuće srednje vrijednosti ulaznih podataka.

Implementacija Metoda najmanjih kvadrata podrazumjeva izračunavanje koeficijenata $a_0$ i $a_1$ iz ulaznih podataka.

Implementacija modela u Pythonu

In [1]:
# Ucitavanje podataka
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (12.0, 9.0)
In [2]:
# Pripremanje i prikaz ulaznih podataka
data = pd.read_csv('data.csv')
X = data.iloc[:, 0]
Y = data.iloc[:, 1]
plt.scatter(X, Y)
plt.show()
In [3]:
# Izračunavanje koeficijenata --> MODEL
X_mean = np.mean(X)
Y_mean = np.mean(Y)

num = 0
den = 0
for i in range(len(X)):
    num += (X[i] - X_mean)*(Y[i] - Y_mean)
    den += (X[i] - X_mean)**2
m = num / den
c = Y_mean - m*X_mean

print ('Koeficijenti su jednaki:')
print('$a_0$= ', c)
print('$a_1$= ',m)
Koeficijenti su jednaki:
('$a_0$= ', 9.908606190326509)
('$a_1$= ', 1.287357370010931)
In [4]:
# Konstrukcija modela --> izracunavnaje Y(x)
Y_model = m*X + c

plt.scatter(X, Y,label='podaci') # stvarne vrijednosti
plt.plot([min(X), max(X)], [min(Y_model), max(Y_model)], label='model',color='red') # model
plt.legend()
plt.show()
In [5]:
## Primjer 2: Matematičko klatno
In [6]:
# Pripremanje i prikaz ulaznih podataka
data = pd.read_csv('data.csv')
X = data.iloc[:, 0]
Y = data.iloc[:, 1]
plt.scatter(X, Y)
plt.show()
In [7]:
# Izračunavanje koeficijenata --> MODEL
X_mean = np.mean(X)
Y_mean = np.mean(Y)

num = 0
den = 0
for i in range(len(X)):
    num += (X[i] - X_mean)*(Y[i] - Y_mean)
    den += (X[i] - X_mean)**2
m = num / den
c = Y_mean - m*X_mean

print ('Koeficijenti su jednaki:')
print('$a_0$= ', c)
print('$a_1$= ',m)
Koeficijenti su jednaki:
('$a_0$= ', 9.908606190326509)
('$a_1$= ', 1.287357370010931)
In [8]:
# Konstrukcija modela --> izracunavnaje Y(x)
Y_model = m*X + c
plt.style.use('Solarize_Light2')
plt.scatter(X, Y,label='podaci') # stvarne vrijednosti
plt.plot([min(X), max(X)], [min(Y_model), max(Y_model)], label='model',color='red') # model
plt.xlabel("l(cm)")
plt.ylabel('T^2(s^2)')
plt.title('Metod najmanjih kvadrata')
plt.legend()
plt.show()

Modeliranje sa SciPy paketom

In [10]:
import scipy as sc
In [11]:
p1=sc.polyfit(X,Y,1)
In [12]:
Y_model2=sc.polyval(p1,X)
In [13]:
plt.style.use('seaborn-darkgrid')
plt.scatter(X, Y,label='podaci') # stvarne vrijednosti
plt.plot([min(X), max(X)], [min(Y_model), max(Y_model)], label='model',color='red') # model
plt.plot([min(X), max(X)], [min(Y_model2), max(Y_model2)],'g-.' ,label='model2') # model2
plt.xlabel("l(cm)")
plt.ylabel('T^2(s^2)')
plt.title('Matematicko klatno')
plt.legend()
plt.show()