Machine Learning 10: Locally weighted Linear Regression - Εκκρεμές
1. Θεωρία: Locally weighted Linear Regression
Η κλασική γραμμική παλινδρόμηση υποθέτει ότι υπάρχει μια ενιαία γραμμική σχέση μεταξύ των μεταβλητών σε όλο το σύνολο των δεδομένων. Δηλαδή το ίδιο μοντέλο χρησιμοποιείται για όλα τα σημεία.
Όταν όμως η πραγματική σχέση μεταξύ X και Y είναι μη γραμμική, ένα μόνο ευθύγραμμο μοντέλο δεν μπορεί να περιγράψει καλά τα δεδομένα.
Η Locally Weighted Linear Regression (LWLR) αντιμετωπίζει αυτό το πρόβλημα κατασκευάζοντας πολλά τοπικά γραμμικά μοντέλα.
Για κάθε σημείο πρόβλεψης, το μοντέλο δίνει μεγαλύτερο βάρος στα κοντινά δεδομένα και μικρότερο βάρος στα πιο μακρινά. Έτσι το μοντέλο προσαρμόζεται τοπικά στη μορφή των δεδομένων.
Με αυτόν τον τρόπο, παρότι κάθε επιμέρους μοντέλο είναι γραμμικό, το συνολικό αποτέλεσμα μπορεί να προσεγγίσει πολύπλοκες μη γραμμικές σχέσεις.
Για αυτό χρησιμοποιούμε βάρη (weights). Κάθε παρατήρηση λαμβάνει ένα βάρος wᵢ.
Η μέθοδος ελαχιστοποιεί το cost function:
Σ wᵢ (yᵢ − (a + b xᵢ))²
Έτσι:
- μεγαλύτερο βάρος → μεγαλύτερη εμπιστοσύνη στη μέτρηση
- μικρότερο βάρος → μικρότερη επίδραση στο αποτέλεσμα
2. Θεωρία Εκκρεμούς
Ένα απλό εκκρεμές αποτελείται από μια μάζα που ταλαντώνεται στο άκρο ενός νήματος μήκους L υπό την επίδραση της βαρύτητας.
Για μικρές γωνίες εκτροπής η περίοδος ταλάντωσης δίνεται από:
T = 2π √(L / g)
- T = περίοδος ταλάντωσης
- L = μήκος νήματος
- g = επιτάχυνση βαρύτητας
Στην ιδανική περίπτωση η περίοδος δεν εξαρτάται από τη γωνία εκτροπής. Όμως για μεγάλες γωνίες εμφανίζεται μη γραμμικότητα και η περίοδος αυξάνεται.
Η πραγματική σχέση προσεγγίζεται από:
T ≈ T₀ (1 + θ²/16)
όπου θ είναι η αρχική γωνία εκτροπής.
Αυτό σημαίνει ότι η σχέση μεταξύ γωνίας και περιόδου δεν είναι γραμμική.
Σε αυτή την περίπτωση ένα απλό γραμμικό μοντέλο δεν είναι αρκετό και χρησιμοποιούμε Locally Weighted Regression.
3. Dataset
Το παρακάτω dataset προσομοιώνει ένα πείραμα εκτροπής εκκρεμούς. Μετράμε την περίοδο ταλάντωσης για διαφορετικές αρχικές γωνίες εκτροπής. Οι μετρήσεις περιέχουν μικρό πειραματικό θόρυβο (measurement noise), όπως συμβαίνει σε πραγματικά εργαστήρια φυσικής.
| ID | Angle (degrees) | Period (seconds) |
|---|---|---|
| 1 | 0 | 2.01 |
| 2 | 2 | 2.00 |
| 3 | 4 | 2.02 |
| 4 | 6 | 2.03 |
| 5 | 8 | 2.01 |
| 6 | 10 | 2.04 |
| 7 | 12 | 2.03 |
| 8 | 14 | 2.05 |
| 9 | 16 | 2.04 |
| 10 | 18 | 2.06 |
| 11 | 20 | 2.07 |
| 12 | 22 | 2.08 |
| 13 | 24 | 2.09 |
| 14 | 26 | 2.11 |
| 15 | 28 | 2.10 |
| 16 | 30 | 2.13 |
| 17 | 32 | 2.12 |
| 18 | 34 | 2.15 |
| 19 | 36 | 2.16 |
| 20 | 38 | 2.18 |
| 21 | 40 | 2.19 |
| 22 | 42 | 2.21 |
| 23 | 44 | 2.23 |
| 24 | 46 | 2.25 |
| 25 | 48 | 2.27 |
| 26 | 50 | 2.29 |
| 27 | 52 | 2.32 |
| 28 | 54 | 2.34 |
| 29 | 56 | 2.37 |
| 30 | 58 | 2.40 |
4. Εκτίμηση μήκους εκκρεμούς και θεωρητικός έλεγχος
Από το dataset παρατηρούμε ότι για μικρές γωνίες εκτροπής η περίοδος είναι περίπου:
T₀ ≈ 2.0 seconds
Για ένα απλό εκκρεμές η περίοδος για μικρές γωνίες δίνεται από:
T = 2π √(L/g)
όπου:
- T = περίοδος ταλάντωσης
- L = μήκος εκκρεμούς
- g = επιτάχυνση βαρύτητας (9.81 m/s²)
Λύνοντας ως προς L:
L = g (T / 2π)²
Αν χρησιμοποιήσουμε την περίοδο για μικρή γωνία (από το dataset):
T ≈ 2.0 s
τότε:
L = 9.81 * (2 / (2π))² L ≈ 0.99 m
Άρα το dataset αντιστοιχεί περίπου σε ένα εκκρεμές μήκους 1 μέτρου.
Θεωρητικός έλεγχος για μεγαλύτερες γωνίες
Για μεγαλύτερες γωνίες εκτροπής η περίοδος δεν είναι πλέον σταθερή και αυξάνεται ελαφρά. Η προσέγγιση για την περίοδο είναι:
T(θ) ≈ T₀ (1 + θ² / 16)
όπου θ είναι η αρχική γωνία εκτροπής σε radians.
Για παράδειγμα, για τη μεγαλύτερη γωνία του dataset:
θ = 58°
Μετατροπή σε radians:
θ = 58 * π / 180 θ ≈ 1.01 rad
Υπολογίζουμε:
T(θ) = 2.0 (1 + θ² / 16) T(58°) = 2.0 (1 + 1.01² / 16) T(58°) ≈ 2.13 s
5. Αλγόριθμος Locally Weighted Regression
Ο παρακάτω κώδικας εφαρμόζει Locally Weighted Linear Regression στο dataset του εκκρεμούς σε python. Τον εκτέλεσα στο Spyder.
import numpy as np
import matplotlib.pyplot as plt
# dataset
angle = np.array([
0,2,4,6,8,10,12,14,16,18,
20,22,24,26,28,30,32,34,36,38,
40,42,44,46,48,50,52,54,56,58
])
period = np.array([
2.01,2.00,2.02,2.03,2.01,2.04,2.03,2.05,2.04,2.06,
2.07,2.08,2.09,2.11,2.10,2.13,2.12,2.15,2.16,2.18,
2.19,2.21,2.23,2.25,2.27,2.29,2.32,2.34,2.37,2.40
])
# design matrix
X = np.vstack([np.ones(len(angle)), angle]).T
y = period
# bandwidth parameter
tau = 10
def lwlr(x_query, X, y, tau):
m = X.shape[0]
W = np.zeros((m, m))
for i in range(m):
diff = x_query - X[i]
W[i, i] = np.exp(-(diff @ diff) / (2 * tau**2))
theta = np.linalg.pinv(X.T @ W @ X) @ (X.T @ W @ y)
return x_query @ theta
# prediction points (smooth curve)
x_plot = np.linspace(0, 58, 100)
X_plot = np.vstack([np.ones(len(x_plot)), x_plot]).T
predictions = np.array([lwlr(x, X, y, tau) for x in X_plot])
# Mean Squared Error on training data
train_pred = np.array([lwlr(x, X, y, tau) for x in X])
mse = np.mean((y - train_pred)**2)
print("MSE:", mse)
# estimate small-angle period
T0 = predictions[0]
print("Estimated small angle period:", T0)
# estimate pendulum length
g = 9.81
L = g * (T0 / (2*np.pi))**2
print("Estimated pendulum length:", L)
# theoretical curve for large angles
theta_rad = np.radians(x_plot)
T_theory = T0 * (1 + theta_rad**2 / 16)
# plot
plt.scatter(angle, period, label="Experimental data")
plt.plot(x_plot, predictions,
linewidth=2,
label="LWLR fit")
plt.xlabel("Angle (degrees)")
plt.ylabel("Period (seconds)")
plt.title("Locally Weighted Regression-Pendulum")
plt.legend()
plt.grid(True)
plt.show()
theta_deg = 58
# ML prediction at 58 degrees
x_query = np.array([1, theta_deg])
T_ml = lwlr(x_query, X, y, tau)
# theoretical prediction
theta_rad = np.radians(theta_deg)
T_theory = T0 * (1 + theta_rad**2 / 16)
# experimental measurement (last dataset point)
T_exp = period[-1]
print("--------------------------")
print("Validation at 58 degrees")
print("--------------------------")
print("Experimental period:", T_exp)
print("ML prediction:", T_ml)
print("Theoretical prediction:", T_theory)
# differences
print()
print("Difference ML vs Experiment:", abs(T_ml - T_exp))
print("Difference Theory vs Experiment:", abs(T_theory - T_exp))
print("Difference ML vs Theory:", abs(T_ml - T_theory))
6. Πώς λειτουργεί το Locally Weighted Linear Regression
Ο αλγόριθμος Locally Weighted Linear Regression (LWLR) δεν υπολογίζει ένα μόνο μοντέλο για όλα τα δεδομένα. Αντίθετα, για κάθε σημείο όπου θέλουμε να κάνουμε πρόβλεψη κατασκευάζει ένα τοπικό γραμμικό μοντέλο.
Η βασική ιδέα είναι ότι τα δεδομένα που βρίσκονται κοντά στο σημείο πρόβλεψης πρέπει να έχουν μεγαλύτερη επίδραση στο αποτέλεσμα σε σχέση με δεδομένα που βρίσκονται μακριά.
Βήμα 1: Δημιουργία του πίνακα δεδομένων
Ο πίνακας σχεδιασμού (design matrix) περιλαμβάνει μια στήλη μονάδων για τον σταθερό όρο του μοντέλου και μια στήλη με τις γωνίες.
X = [ 1 angle₁
1 angle₂
1 angle₃
...
1 angleₙ ]
Το μοντέλο που θέλουμε να προσαρμόσουμε είναι:
T = a + b·θ
όπου:
- a = σταθερός όρος (intercept)
- b = κλίση της ευθείας
- θ = γωνία εκτροπής
Βήμα 2: Υπολογισμός των βαρών (weights)
Για κάθε σημείο πρόβλεψης υπολογίζουμε ένα βάρος για κάθε παρατήρηση του dataset. Τα βάρη υπολογίζονται με μια συνάρτηση Gaussian kernel.
wᵢ = exp( − ||x_query − xᵢ||² / (2τ²) )
όπου:
- x_query = το σημείο όπου θέλουμε πρόβλεψη
- xᵢ = το i-οστό δεδομένο του dataset
- τ = παράμετρος bandwidth
Η παράμετρος τ (tau) καθορίζει πόσο τοπικό είναι το μοντέλο.
- μικρό τ → επηρεάζουν κυρίως τα πολύ κοντινά σημεία
- μεγάλο τ → επηρεάζουν περισσότερα σημεία
Στον κώδικα χρησιμοποιούμε:
tau = 10
Τα βάρη αποθηκεύονται σε έναν διαγώνιο πίνακα:
W = [w₁ 0 0 ... 0 0 w₂ 0 ... 0 0 0 w₃ ... 0 ... 0 0 0 ... wₙ]
Βήμα 3: Υπολογισμός των παραμέτρων του τοπικού μοντέλου
Για να βρούμε τις παραμέτρους του τοπικού γραμμικού μοντέλου χρησιμοποιούμε την παρακάτω εξίσωση (weighted least squares):
θ = (Xᵀ W X)⁻¹ Xᵀ W y
όπου:
- X = πίνακας δεδομένων
- W = πίνακας βαρών
- y = μετρήσεις της περιόδου
Στον κώδικα αυτό υλοποιείται ως (στην Python το @ χρησιμοποιείται στον πολλαπλασιασμό πινάκων):
theta = np.linalg.pinv(X.T @ W @ X) @ (X.T @ W @ y)
Χρησιμοποιούμε την ψευδοαντίστροφη (pseudo-inverse) ώστε ο υπολογισμός να είναι αριθμητικά σταθερός.
Βήμα 4: Υπολογισμός της πρόβλεψης
Αφού υπολογίσουμε τις παραμέτρους του μοντέλου, η πρόβλεψη για το σημείο x_query είναι:
prediction = x_query · θ
Η διαδικασία επαναλαμβάνεται για πολλά σημεία, ώστε να δημιουργηθεί μια ομαλή καμπύλη.
Γιατί το μοντέλο γίνεται μη γραμμικό
Κάθε πρόβλεψη χρησιμοποιεί διαφορετικά βάρη και επομένως διαφορετικό τοπικό γραμμικό μοντέλο.
Έτσι, παρόλο που κάθε επιμέρους μοντέλο είναι γραμμικό, η συνολική καμπύλη μπορεί να προσεγγίσει πολύπλοκες μη γραμμικές σχέσεις.
7. Αποτελέσματα
MSE: 0.00010464897425655976 Estimated small angle period: 2.0024417220910875 Estimated pendulum length: 0.996389269102797 -------------------------- Validation at 58 degrees -------------------------- Experimental period: 2.4 ML prediction: 2.3883843777333515 Theoretical prediction: 2.130689729186646 Difference ML vs Experiment: 0.011615622266648362 Difference Theory vs Experiment: 0.2693102708133539 Difference ML vs Theory: 0.25769464854670554
Το σφάλμα Mean Squared Error (MSE) είναι περίπου 0.00010, τιμή που δείχνει ότι η καμπύλη που υπολογίζει το μοντέλο βρίσκεται πολύ κοντά στις πειραματικές μετρήσεις.
Από το regression εκτιμήθηκε περίοδος μικρής γωνίας T₀ ≈ 2.002 s, η οποία είναι πολύ κοντά στη θεωρητική τιμή των 2 s. Χρησιμοποιώντας αυτή την τιμή υπολογίστηκε μήκος εκκρεμούς περίπου L ≈ 0.996 m, δηλαδή πολύ κοντά στο πραγματικό μήκος του 1 m.
Για τη μεγάλη γωνία 58° το μοντέλο προβλέπει περίοδο 2.388 s, ενώ η πειραματική μέτρηση είναι 2.40 s. Η διαφορά είναι μόλις 0.012 s, γεγονός που δείχνει ότι το μοντέλο προσαρμόζεται πολύ καλά στα δεδομένα του πειράματος.
Βιβλιογραφία
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2021). Introduction to linear regression analysis (6th ed.). Wiley.
Stanford Online. Stanford CS229: Machine Learning - Locally Weighted & Logistic Regression | Lecture 3 (Autumn 2018) [Video]. YouTube. https://youtu.be/het9HFqo1TQ?list=PLoROMvodv4rMiGQp3WXShtMGgzqpfVfbU
St-Aubin, A. (n.d.). An introduction to supervised machine learning. McGill University.