Biblioteka numpy zawiera procedury do aproksymacji. W szczególności:
numpy.linalg.lstsq
,numpy.polyfit
,import matplotlib.pyplot as plt
import numpy as np
#from scipy.interpolate import interp1d
#from scipy.interpolate import lagrange
#from scipy.interpolate import BarycentricInterpolator
Plik AVERAGE300aprox.dat
zawiera dane już oczyszczone.
data = np.loadtxt("AVERAGE300aprox.dat")
plt.plot(data[:,0],data[:,1],'.')
plt.show()
Dokonamy najprostszej aproksymacji wielomianem stopnia trzeciego, używając funkcji polyfit()
.
Funkcja zwraca współczynniki wielomianu (w przypadku gdy stopień, jak w przykładzie jest równy 3) $$z_0 x^3 + z_1 x^2 + z_2 x + z_3$$
z = np.polyfit(data[:,0],data[:,1], 3)
Jak widać funkcja informuje, że zadanie jest źle uwarunkowane. Wynika to (w tym przypadku) z tego, że wartości $x$ są bardzo duże (gdyż reprezentują czas bezwzględny w sekundach). Bardzo łatwo zamienić je na wartości reprezentujące czas względny (względem początku doby) w sekundach odejmując od każdej wartości odciętej minimum.
tmin = min(data[:,0])
data[:,0] = data[:,0] - tmin
z = np.polyfit(data[:,0],data[:,1], 3)
Jak widać jest wyraźnie lepiej (nie pojawia się komunikat).
z
jest tablicą zwierającą współczynniki wielomianu. Zwracam uwagę, że ze względu na spore wartości $x$ (liczba sekund w ciągu doby) współczynniki są bardzo małe. Grozi to utratą dokładności podczas obliczeń. Można spróbować zamienić sekundy na godziny dzieląc odcięte przez 3600.
Jeszcze lepsze efekty osiągniemy normalizując odcięte danych do zakresu $[-1, 1]$. Daje to szansę, że skrajne wartości podnoszone do dużych nawet potęg będą przyjmowały sensowniejsze wartości.
z
data[:,0] = data[:,0]/3600.
z = np.polyfit(data[:,0],data[:,1], 3)
z
Współczynniki wielomianu uległy modyfikacji, za wyjątkiem ostatniego $z_3$ reprezentującego wartość funkcji w punkcie 0.
Do manipulacji na wielomianie warto używac funkcji poly1d()
wielomian = np.poly1d(z)
wielomian(0.)
wielomian(24)
Obejrzyjmy na wykresie rezultat aproksymacji.
x=np.linspace(0, 24, 25)
plt.plot(x, wielomian(x),'g-',data[:,0], data[:,1], 'r.')
plt.show()
Przybliżenie nie jest idealne, ale (w pewnym sense) oddaje charakter dobowych zmian temperatury: w ciągu doby rosną, a na zakończenie dnia temperatura jest mniejsza niż na jego początku.
Niech N
będzie to stopień wielomianu interpolującego. W poniższym przykłądzie można dowlonie(?) zmieniać N
i przeliczać…
N = 18
plt.plot(x, np.poly1d(np.polyfit(data[:,0],data[:,1], N))(x),'g-',
data[:,0], data[:,1], 'r.')
plt.show()
Zobaczmy jak wyglądają różnice wartości prawdziwej i aproksymowanej.
err = data[:,1] - np.poly1d(np.polyfit(data[:,0],data[:,1], N))(data[:,0])
plt.plot(data[:,0], err)
plt.show()
Generalnie nie jest źle: wartości skupiają się wokół zera w sposób dosyć symetryczny. Zobaczmy histogram.
plt.hist(err)
plt.show()
Nie do końca przypomina on histogram rozkłądu normalnego, ale nie jest ,,najgorszy''.
Można też zdefiniować sobie funkcję (na przykład wykres
) i korzystać z niej w celu zobaczenia rezultatu aproksymacji.
def wykres(N):
plt.plot(x, np.poly1d(np.polyfit(data[:,0],data[:,1], N))(x),'g-',
data[:,0], data[:,1], 'r.')
plt.show()
wykres(5)
Można też zabawić się w wykres interakcyjny, w którym suwaczkiem zmieniamy wartość parametru. Ale to już jakaś magia…
%matplotlib inline
#from __future__ import print_function
from ipywidgets import interactive
import matplotlib.pyplot as plt
import numpy as np
def wykres_wykres(N):
plt.figure(2)
plt.plot(x, np.poly1d(np.polyfit(data[:,0],data[:,1], N))(x),'g-',
data[:,0], data[:,1], 'r.')
# plt.ylim(-5, 5)
plt.show()
interactive_plot = interactive(wykres_wykres, N=(0, 30))
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot
Dla $N=0$ dostajemy średnią temperaturę dobową, a dla $N=1$ trend — temperatura spada.
Również dla kolejnych wartości $N$ można znaleźć jakąś interpretację…
Niestety, dla $N >18$ funkcja informuje, że zadanie staje się źle uwarunkowane. Jest to znowu związane z tym, że wartości $x$ podniesiona do potęgi 19 czy 20 staje się baaaardzo duża:
24.**20
Spójrzmy zatem na to co składa się na wartość wielomainu dla czasu $x$.
N = 19
x = 24
a=np.polyfit(data[:,0],data[:,1], N)
for i in range(0, N):
print(a[i]*x**(N-i))
Jak widać sumowane są liczby o bardzo dużej rozpiętości wartości. Musi prowadzić to do problemów…