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)
/home/myszka/anaconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3418: RankWarning: Polyfit may be poorly conditioned exec(code_obj, self.user_global_ns, self.user_ns)
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
array([-1.50077850e-15, -4.89454863e-09, 3.67306502e-04, -2.29165068e+00])
data[:,0] = data[:,0]/3600.
z = np.polyfit(data[:,0],data[:,1], 3)
z
array([-7.00203219e-05, -6.34333503e-02, 1.32230341e+00, -2.29165068e+00])
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.)
-2.2916506803099526
wielomian(24)
-8.06193960023452
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
interactive(children=(IntSlider(value=15, description='N', max=30), Output(layout=Layout(height='350px'))), _d…
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
4.0199887178406037e+27
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)
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:3: RankWarning: Polyfit may be poorly conditioned This is separate from the ipykernel package so we can avoid doing imports until
for i in range(0, N):
print(a[i]*x**(N-i))
-11330548833.759817 91343839406.0988 -323839414937.3468 649910659240.3022 -768190706710.5015 429110280561.0846 186788204073.06516 -623109225162.7025 652222465812.5778 -426061915498.8431 193788821235.2673 -63305881981.175026 14896926263.546213 -2489387080.6460786 286464930.539979 -21521206.64991253 960976.6069019116 -21266.746996815928 173.36912084905677
Jak widać sumowane są liczby o bardzo dużej rozpiętości wartości. Musi prowadzić to do problemów…