5. Παρεμβολή και παλινδρόμηση

5.1. Διατύπωση του προβλήματος

Πολλές φορές αντί της μαθηματικής συνάρτησης \(y=f(x)\) μεταξύ δύο μεγεθών \(x\) και \(y\) είναι γνωστή μία σειρά από \(n\) σημεία \((x_i,y_i)\), η οποία περιγράφει προσεγγιστικά την εξάρτηση. Χαρακτηριστικό παράδειγμα είναι οι πίνακες που βρίσκονται στα παραρτήματα επιστημονικών βιβλίων και τεχνικών εγχειριδίων, καθώς και τα αποτελέσματα προσομοιώσεων. Η συνάρτηση \(f(x)\) μπορεί να μας είναι άγνωστη, να μην υπάρχει αναλυτική λύση της ή να είναι πολύ δύσκολος ο υπολογισμός της. Για την εύρεση ενός \(y\) σε τυχαία θέση \(x\) μπορεί να χρησιμοποιηθεί η σειρά των σημείων χρησιμοποιώντας μια μέθοδο παρεμβολής.

Exercise 5.1

Μια μετεωρολογική υπηρεσία ανακοινώνει τις προγνώσεις της με χρονικό βήμα 3 ωρών. Στον Πίνακα 5.1 φαίνεται η πρόβλεψη της θερμοκρασίας για ένα εικοσιτετράωρο. Απεικονίστε την θερμοκρασία σε γράφημα με σημεία (χωρίς γραμμές) και εκτιμήστε την θερμοκρασία στις 05:00 και 20:00.

Table 5.1 Πρόβλεψη θερμοκρασίας.

Πολιτική ώρα [h]	Θερμοκρασία [°C]
0	18
3	16
6	15
9	21
12	25
15	27
18	26
21	23
24	18

Solution to Exercise 5.1

import numpy as np
import matplotlib.pyplot as plt

x = np.array([0, 3, 6, 9, 12, 15, 18, 21, 24])
y = np.array([18, 16, 15, 21, 25, 27, 26, 23, 18])

fig, ax = plt.subplots()
ax.set(xlabel="Time [h]", ylabel="Temperature [°C]")

ax.plot(x, y, "o")
plt.xlim(0, 24)
plt.ylim(12, 28)
ax.grid()
plt.show()

Από το γράφημα φαίνεται ότι στις 05:00 η θερμοκρασία θα είναι περίπου 15.5°C και στις 20:00 24°C.

5.2. Γραμμική παρεμβολή

Η πιο απλή και συνηθισμένη μέθοδος παρεμβολής είναι η γραμμική παρεμβολή κατά την οποία θεωρούμε ότι ευθύγραμμα τμήματα ενώνουν τα γνωστά μας σημεία. Το ευθύγραμμο τμήμα σε κάθε διάστημα εκφράζεται μαθηματικά ως:

\[ y=y_i+\frac{y_{i+1}-y_i}{x_{i+1}-x_i}(x-x_i) \;\text{για}\;x_i<x<x_{i+1}\]

Exercise 5.2

Λύστε την Άσκηση 5.1 εφαρμόζοντας γραμμική παρεμβολή με αναλυτικές σχέσεις.

Solution to Exercise 5.2

Πρώτα θα προσδιορίσουμε τα σημεία, τα οποία θα χρησιμοποιήσουμε για τον υπολογισμό. Για τον υπολογισμό του χρόνου \(x=5\) θα χρησιμοποιήσουμε τα κοντινότερα σημεία, δηλαδή αυτά με \(x_i=3\) και \(x_{i+1}=6\), επομένως:

\[ y= 16+\frac{16-15}{6-3}(5-3)=15.33\]

Ομοίως για \(x=20\) στο αντίστοιχο διάστημα.

import numpy as np
import matplotlib.pyplot as plt

x = np.array([0, 3, 6, 9, 12, 15, 18, 21, 24])
y = np.array([18, 16, 15, 21, 25, 27, 26, 23, 18])

fig, ax = plt.subplots()
ax.set(xlabel="Time [h]", ylabel="Temperature [°C]")

i = 1
xa = 5
ya = y[i] + (y[i + 1] - y[i]) / (x[i + 1] - x[i]) * (xa - x[i])
print(xa, ya)

i = 6
xb = 20
yb = y[i] + (y[i + 1] - y[i]) / (x[i + 1] - x[i]) * (xb - x[i])
print(xb, yb)

ax.plot(x, y, "o-")
ax.plot([xa, xb], [ya, yb], "rD")  # red diamonds
plt.xlim(0, 24)
plt.ylim(12, 28)
ax.grid()
plt.show()

5 15.333333333333334
20 24.0

Η γραμμική παρεμβολή δίνει καλά αποτελέσματα όταν η παράσταση ενός μεγέθους προσεγγίζεται ικανοποιητικά από ευθύγραμμα τμήματα. Σε περιπτώσεις με μεγάλη καμπυλότητα απαιτούνται πιο πολλά σημεία για ικανοποιητική ακρίβεια ή μια καλύτερη συνάρτηση παρεμβολής. Η γραμμική παρεμβολή δύναται να εφαρμοστεί και εκτός του διαστήματος των γνωστών σημείων \(x_i\) (extrapolation). Η προσέγγιση αυτή συναντάται με την ονομασία προεκβολή, παρεκβολή ή και παρέκταση. Η προεκβολή πρέπει να χρησιμοποιείται με προσοχή, καθώς η ακρίβειά της μειώνεται μακρυά από τα γνωστά σημεία.

Important

Η συνάρτηση numpy.interp δέχεται τουλάχιστον τρία ορίσματα, το σημείο της παρεμβολής \(x\) και δύο διανύσματα ίδιου μεγέθους με τα γνωστά σημεία \(x_i\) και \(y_i\). Το διάνυσμα των \(x_i\) πρέπει να έχει αύξουσα σειρά, έτσι ώστε η συνάρτηση να εντοπίσει το κατάλληλο διάστημα \(x_i,x_{i+1}\) που περικλείει το \(x\) της παρεμβολής. Από προεπιλογή η συνάρτηση δεν κάνει προεκβολή, αλλά επιστρέφει την τιμή της \(y_1\) αν \(x<x_1\) και την \(y_n\) αν \(x>x_n\). Το \(x\) μπορεί να είναι κι αυτό διάνυσμα, οπότε κι η συνάρτηση θα επιστρέψει ένα διάνυσμα με τιμές \(y\).

Exercise 5.3

Λύστε την Άσκηση 5.1 εφαρμόζοντας γραμμική παρεμβολή με την συνάρτηση numpy.interp.

Solution to Exercise 5.3

import numpy as np

x = np.array([0, 3, 6, 9, 12, 15, 18, 21, 24])
y = np.array([18, 16, 15, 21, 25, 27, 26, 23, 18])

xa = 5
ya = np.interp(xa, x, y)
print(xa, ya)

xb = 20
yb = np.interp(xb, x, y)
print(xb, yb)

xc = 27
yc = np.interp(xc, x, y)  # no extrapolation
print(xc, yc)

5 15.333333333333334
20 24.0
27 18.0

5.3. Παρεμβολή πολυωνύμου

Στο τελευταίο γράφημα φαίνεται ότι τα ευθύγραμμα τμήματα σχηματίζουν γωνίες μεταξύ τους, οι οποίες φαίνονται πιο καθαρά στις περιοχές που αναμένεται μεγαλύτερη καμπυλότητα. Τα ευθύγραμμα τμήματα δίνουν ασυνεχή πρώτη παράγωγο (κλίση) στα εσωτερικά σημεία, γεγονός μη αναμενόμενο. Για να έχουμε συνέχεια μεταξύ των τμημάτων μπορούμε να παρεμβάλουμε ένα πολυώνυμο που περνάει από όλα τα σημεία. Το πολυώνυμο αυτό θα πρέπει να είναι το λιγότερο \(n-1\) βαθμού:

\[f(x)=a_0+a_1x+a_2x^2+ \ldots +a_{n-1}x^{n-1} \]

Ζητούμενο είναι να προσδιοριστούν οι \(n\) συντελεστές \(a_{i}\) του πολυωνύμου. Από τα \(n\) γνωστά σημεία προκύπτει ένα γραμμικό σύστημα \(n\) εξισώσεων κατάλληλο για την εύρεση των άγνωστων συντελεστών. Το σύστημα αυτό μπορεί να επιλυθεί με τις μεθόδους που αναφέρθηκαν στο Κεφάλαιο 4. Μπορούν όμως να χρησιμοποιηθούν και πολυώνυμα ειδικής μορφής που ελαχιστοποιούν τις μαθηματικές πράξεις.

5.3.1. Πολυώνυμα Newton

Τα πολυώνυμα Newton έχουν την μορφή:

\[ f(x) = a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + \dots + a_{n-1}(x-x_0)(x-x_1)\dots(x-x_{n-1}) \]

ή ισοδύναμα:

\[ f(x) = \sum_{i=0}^{n-1}{a_in_i(x)} \]

όπου:

\[n_i(x) = \prod_{j=0}^{i}(x-x_j)\]

Αν και η μορφή αυτή είναι πιο πολύπλοκη μορφή από την γενική που παρουσιάστηκε παραπάνω, έχει το πλεονέκτημα του εύκολου υπολογισμού των συντελεστών \(a_i\). Στον υπολογισμό χρησιμοποιούνται οι πεπερασμένες διαφορές που ορίζονται ως:

\[\begin{split} \begin{align} f[x_i]&=&y_i\\ f[x_{i+1},x_{i}]&=&\frac{f[x_{i+1}]-f[x_i]}{x_{i+1}-x_{i}}&=&\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\\ f[x_i,x_{i-1},\dots,x_1,x_0]&=&\frac{f[x_i,\dots,x_1]-f[x_{i-1},x_0]}{x_i-x_0} \end{align} \end{split}\]

Οι πεπερασμένες διαφορές ανώτερης τάξης υπολογίζονται με αυτές της αμέσως προηγούμενης τάξης, επομένως μπορεί να εφαρμοστεί η παρακάτω ακολουθία στήλη προς στήλη:

\[\begin{split} \begin{align} x_0 && \pmb{f[x_0]} \\ && &\rangle & \pmb{f[x_1,x_0]} \\ x_1 && f[x_1] & & & \rangle & \pmb{f[x_2, x_1,x_0]}\\ && &\rangle &f[x_2,x_1] & & & \rangle & \pmb{f[x_3, x_2, x_1,x_0]}\\ x_2 && f[x_2] & & & \rangle & f[x_3, x_2,x_1] & & &\rangle & \pmb{f[x_4, x_3, x_2, x_1,x_0]}\\ && &\rangle &f[x_3,x_2] & & & \rangle & f[x_4, x_3, x_2, x_1]\\ x_3 && f[x_3] & & & \rangle & f[x_4, x_3,x_2]\\ && &\rangle & f[x_4,x_3] \\ x_4 && f[x_4] \end{align} \end{split}\]

Οι συντελεστές του πολυωνύμου Newton είναι οι κορυφαίες τιμές κάθε στήλης (σημειώνονται με έντονη γραφή), δηλαδή:

\[\begin{split} a_0=f[x_0]\\ a_1=f[x_1,x_0]\\ \dots \end{split}\]

Προσθέτοντας ένα ακόμα σημείο (το \(x_5\) στο παράδειγμα) υπολογίζεται μια γραμμή στην βάση και επεκτείνεται η διάταξη κατά μία στήλη. Έτσι μπορεί εύκολα να μεγαλώσει ο βαθμός του πολυωνύμου παρεμβολής, ένα χρήσιμο χαρακτηριστικό όταν δεν γνωρίζουμε εκ των προτέρων τον επιθυμητό βαθμό.

5.3.2. Πολυώνυμα Lagrange

Μια παρόμοια μέθοδος είναι τα πολυώνυμα Lagrange, η οποία διευκολύνει κι αυτή τους υπολογισμούς. Τα πολυώνυμα Lagrange έχουν την μορφή:

\[L(x)=\sum_{i=1}^{i=n} y_i L_i(x) \]

όπου τα γινόμενα \(L_i\):

\[L_i(x)=\prod_{j=1,j\neq i}^{j=n} \frac{x-x_j}{x_i-x_j} \]

Το \(j\) παίρνει όλες τις τιμές από i έως n εκτός από το \(i\).

Παρατηρείστε ότι:

\[\begin{split} L_i(x_0)=0\\ L_i(x_1)=0\\ \dots\\ L_i(x_i)=1\\ \dots \end{split}\]

H SciPy περιλαμβάνει την σχετική συνάρτηση scipy.interpolate.lagrange.

Exercise 5.4

Λύστε την Άσκηση 5.1 εφαρμόζοντας πολυωνυμική παρεμβολή με την συνάρτηση scipy.interpolate.lagrange.

Solution to Exercise 5.4

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import lagrange

x = np.array([0, 3, 6, 9, 12, 15, 18, 21, 24])
y = np.array([18, 16, 15, 21, 25, 27, 26, 23, 18])

fig, ax = plt.subplots()
ax.set(xlabel="Time [h]", ylabel="Temperature [°C]")

f = lagrange(x, y)
xf = np.linspace(0, 24, 24 * 10)

ax.plot(x, y, "o-")
ax.plot(xf, f(xf), "-")
ax.plot([5, 20], [f(5), f(20)], "rD")  # red Diamonds
plt.xlim(0, 24)
plt.ylim(12, 28)
ax.grid()
plt.show()

print(f(5), f(20))

13.830869278058355 23.834527257281152

Important

Πρέπει να σημειωθεί ότι τα πολυώνυμα ικανοποιούν ακριβώς όλα τα σημεία. Ταυτόχρονα είναι πολύ ευαίσθητα σε σφάλματα που μπορεί να εμπεριέχουν τα δεδομένα μας με συνέπεια να παρουσιάζουν πολλαπλά μέγιστα και ελάχιστα ή γενικώς μη αναμενόμενη συμπεριφορά από φυσικής άποψης. Αυτό συμβαίνει στο διάστημα 0-6 h του παραδείγματος μας. Συνήθως η πιο απλή συμπεριφορά είναι πιο κοντά στην πραγματικότητα.

5.4. Παρεμβολή spline

Για την αποφυγή των μειονεκτημάτων της γραμμικής παρεμβολής και της παρεμβολής πολυωνύμων, αναπτύχθηκαν τα spline. Παρόμοια με την γραμμική παρεμβολή, το πεδίο χωρίζεται σε τμήματα στα οποία χρησιμοποιείται διαφορετική συνάρτηση η παρεμβολής. Επιπλέον η συνάρτηση αυτή είναι συνήθως πολυώνυμο τρίτου βαθμού, έτσι ώστε να διασφαλίζεται η συνέχεια στην πρώτη και δεύτερη παράγωγο.
Στη μονάδα κώδικα SciPy ορίζεται ένα τμήμα για κάθε διάστημα μεταξύ δύο σημείων. Η ίδια προσέγγιση χρησιμοποιείται εδώ, χωρίς να είναι αυτή η μόνη επιλογή.

Επομένως για τα \(n-1\) διαστήματα πρέπει να οριστούν \(n-1\) εξισώσεις της μορφής:

\[ S_i(x) = a_i x^3 + b_i x^2 + c_i x + d_i\;\text{για}\;x_{i-1}<x<x_{i} \]

Κάθε πολυώνυμο έχει 4 άγνωστους συντελεστές άρα πρέπει να καταστρωθεί ένα σύστημα \(4(n-1)\) εξισώσεων:

1. Οι εξισώσεις παρεμβολής πρέπει να δίνουν τις ήδη γνωστές τιμές στα άκρα τους

\[\begin{split} \begin{align} S_i(x_i)&=&y_i& \;\text{για}\;i=0,n-2\;\text{(αριστερό σημείο)}\\ S_i(x_{i+1})&=&y_{i+1}& \;\text{για}\;i=0,n-2\;\text{(δεξιό σημείο)} \end{align} \end{split}\]

\[\bf{2(n-1)}\;\textbf{εξισώσεις}\]

2. Οι συναρτήσεις παρεμβολής από της δύο πλευρές ενός εσωτερικού σημείου πρέπει να έχουν ίσες πρώτες και δεύτερες παραγώγους

\[\begin{split} \begin{align} S'_i(x_{i})&=&S'_{i+1}(x_{i}) &\;\text{για}\;i=1,n-2\;\text{(πρώτη παράγωγος)}\\ S''_i(x_{i})&=&S''_{i+1}(x_{i}) &\;\text{για}\;i=1,n-2\;\text{(δεύτερη παράγωγος)} \end{align} \end{split}\]

\[\bf{2(n-2)}\;\textbf{εξισώσεις}\]

3. Τέλος επιβάλουμε αυθαίρετα δύο ακόμα περιορισμούς για να οριστεί πλήρως το σύστημα. Μια συνηθισμένη προσέγγιση είναι η μηδενική καμπυλότητα στα ακραία σημεία \(x_i\)

\[\begin{split} \begin{align} S''_1(x_{0})&=&0&\;\text{(πρώτο σημείο)}\\ S''_{n-1}(x_{n-1})&=&0&\;\text{(τελευταίο σημείο)} \end{align} \end{split}\]

\[\bf{2}\;\textbf{εξισώσεις}\]

Tip

Η πλήρης ονομασία της παραπάνω προσέγγισης είναι παρεμβολή με φυσικά κυβικά spline. Είναι διαθέσιμη στην συνάρτηση scipy.interpolate.CubicSpline με την επιλογή φυσικών οριακών συνθηκών για τα δύο άκρα (bc_type=’natural’). Εξορισμού bc_type=”not-a-knot”, που οδηγεί στην χρήση κοινής συνάρτησης παρεμβολής στα 2 πρώτα και τα 2 τελευταία τμήματα.

Exercise 5.5

Λύστε την Άσκηση 5.1 χρησιμοποιώντας φυσικά κυβικά spline (της μονάδας κώδικα SciPy) και σχολιάστε τις διαφορές στην μορφή της καμπύλης στα αποτελέσματα.

Solution to Exercise 5.5

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline

x = np.array([0, 3, 6, 9, 12, 15, 18, 21, 24])
y = np.array([18, 16, 15, 21, 25, 27, 26, 23, 18])

fig, ax = plt.subplots()
ax.set(xlabel="Time [h]", ylabel="Temperature [°C]")

f = CubicSpline(x, y, bc_type="natural")
xf = np.linspace(0, 24, 24 * 10)

ax.plot(x, y, "o-")
ax.plot(xf, f(xf), "-")
ax.plot([5, 20], [f(5), f(20)], "rD")  # red Diamonds
plt.xlim(0, 24)
plt.ylim(12, 28)
ax.grid()
plt.show()

print(f(5), f(20))

14.647248131784215 24.230036000654557

Παρατηρούμε ότι η διακύμανση της θερμοκρασίας φαίνεται πιο φυσική με χρήση spline. Οι τιμές στις ζητούμενες χρονικές στιγμές διαφοροποιούνται αισθητά από την γραμμική παρεμβολή. Ειδικά στο χρονικό διάστημα 3-6h η γραμμική συνάρτηση δεν ταιριάζει με ρεαλιστική θερμοκρασιακή συμπεριφορά σε αντίθεση με την παρεμβολή spline.

5.5. Παλινδρόμηση

Όταν τα δεδομένα μας περιέχουν σφάλματα μπορούμε να χρησιμοποιήσουμε μια συνάρτηση \(f(x)\) που δεν περνάει απαραίτητα από τα \(m\) σημεία μας, αλλά αρκετά κοντά. Σε αυτή την περίπτωση δεν εφαρμόζουμε μεθόδους παρεμβολή αλλά παλινδρόμησης (regression) και η επιλογή της κατάλληλης συνάρτησης γίνεται ακόμα πιο σημαντική. Εδώ θα εξετάσουμε την γενική περίπτωση πολυωνύμου \(n\) βαθμού, η οποία περιλαμβάνει και την γραμμική παλινδρόμηση:

\[f(x)=a_0+a_1x+a_2x^2+ \ldots +a_{n}x^{n} \]

Οι \(n+1\) συντελεστές \(a_j\) πρέπει να επιλεγούν έτσι ώστε να ελαχιστοποιούνται οι αποκλίσεις \(E_i=y_i-f(x_i)\) στα γνωστά σημεία. Η πιο συνηθισμένη προσέγγιση είναι να ελαχιστοποιείται το άθροισμα των τετραγώνων των σφαλμάτων, δηλαδή η μαθηματική έκφραση:

\[S_{r}=\sum_{i=1}^{m} [y_i-f(x_i)]^2=\sum_{i=1}^{m} [y_i-a_0-a_1{x_i}-a_2{x_i}^2- \ldots -a_{n}{x_i}^{n}]^2\]

Το μέγεθος αυτό ελαχιστοποιείται εκεί που οι μερικές παράγωγοί του ως προς \(a_j\) μηδενίζονται, δηλαδή:

\[\begin{split} \begin{align} \frac{\partial S_r}{\partial a_0}&=&-2\sum_{i=1}^{n} [y_i-a_0-a_1{x_i}-a_2{x_i}^2- \ldots -a_{n}{x_i}^{n}]1\\ \frac{\partial S_r}{\partial a_1}&=&-2\sum_{i=1}^{n} [y_i-a_0-a_1{x_i}-a_2{x_i}^2- \ldots -a_{n}{x_i}^{n}]{x_i}\\ \dots\\ \frac{\partial S_r}{\partial a_n}&=&-2\sum_{i=1}^{n} [y_i-a_0-a_1{x_i}-a_2{x_i}^2- \ldots -a_{n}{x_i}^{n}]{x_i}^n \end{align} \end{split}\]

Έτσι προκύπτει ένα γραμμικό σύστημα, το οποίο μπορεί να λυθεί εύκολα.

Η NumPy περιέχει μια ειδική κλάση για τον χειρισμό πολυωνύμων. Για την εύρεση των συντελεστών διατίθεται η μέθοδος fit, η οποία δέχεται ως ορίσματα τα διανύσματα \(x\) και \(y\) και τον βαθμό του πολυωνύμου.

Για την αξιολόγηση της προσέγγισης του πολυωνύμου (αλλά και οποιασδήποτε συνάρτησης παρεμβολής) στα δεδομένα μας χρησιμοποιείται συχνά ο συντελεστής προσδιορισμού (\(R^2\)):

\[R^2=1-\frac{S_{r}}{S_{t}}=1-\frac{\sum_{i=1}^{n} [y_i-f(x_i)]^2}{\sum_{i=1}^{n} [y_i-\bar{y}]^2}\]

Όπως φαίνεται από τον κλάσμα, όταν οι αποκλίσεις τείνουν στο μηδέν, το \(R^2\) τείνει στην μονάδα. Αντιθέτως όταν η συνάρτηση μας είναι το ίδιο καλή με μια μέση τιμή \(\bar{y}\), δηλαδή δεν εξηγεί καθόλου την εξάρτηση από το \(x\), το \(R^2\) τείνει στο 0.

Exercise 5.6

Εφαρμόστε πολυωνυμική παλινδρόμηση 5ου βαθμού στα δεδομένα της Άσκησης 5.1 και υπολογίστε την θερμοκρασία στα δύο χρονικά σημεία. Αναπαραστήστε τα αποτελέσματα σε γράφημα στο οποίο αναγράφεται κι ο συντελεστής προσδιορισμού.

Solution to Exercise 5.6

import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial

x = np.array([0, 3, 6, 9, 12, 15, 18, 21, 24])
y = np.array([18, 16, 15, 21, 25, 27, 26, 23, 18])

fig, ax = plt.subplots()
ax.set(xlabel="Time [h]", ylabel="Temperature [°C]")

f = Polynomial.fit(x, y, deg=5)
xf = np.linspace(0, 24, 24 * 10)

Rsq = (
    np.corrcoef(y, f(x))[0, 1] ** 2
)  # Calculation of R2 from Pearson correlation matrix.

print(Rsq)

ax.plot(x, y, "o-")
ax.plot(xf, f(xf), "-")
ax.plot([5, 20], [f(5), f(20)], "rD")  # red diamonds
plt.xlim(0, 24)
plt.ylim(12, 28)
plt.text(20, 27, f"R²={Rsq:0.1%}")

ax.grid()
plt.show()

print(f(5), f(20))

0.9771452505827505

15.568138651471992 24.107744107744125