NumPy-polynomit

Kemian tekniikassa haluamme usein sovittaa polynomeja mittausdataan.  Polynomien käsittelyssä voimme käyttää numpy.poly1d-toimintoa (varsinaisesti kyseessä on "luokka", joita käsitellään olio-ohjelmoinnin yhteydessä kurssin 6. kierroksella).

Polynomien luominen

import numpy as np
import matplotlib.pyplot as plt

# Luodaan polynomi x^2 + 4x + 3
# Kutsutaan poly1d:tä kertoimilla [1, 4, 3], eli 1*x^2 + 4*x^1 + 3*x^0
pol = np.poly1d([1, 4, 3]) 
print("Polynomi on (eksponentit ylärivillä):\n", pol)

tulostaa

Polynomi on (eksponentit ylärivillä):
    2
1 x + 4 x + 3

Polynomien perusominaisuudet

# Polynomin arvon laskeminen yksittäisessä pisteessä
pol = np.poly1d([1, 4, 3])
print("Arvo pisteessa x = 5:", pol(5)) 
# Laskee siis 5^2 + 5*4 + 3 = 48

# Polynomin juuret (nollakohdat): pol.r
print("Juuret:", pol.r)

# Polynomin derivaatta: pol.deriv
# pol.deriv palauttaa uuden polynomin (poly1d-olion )
der = pol.deriv()
print("Derivaatta:", der)
print("Derivaatan arvo pisteessä x = 5:", der(5))

tulostaa

Arvo pisteessa x = 5: 48
Juuret: [-3. -1.]
Derivaatta:  
2 x + 4
Derivaatan arvo pisteessä x = 5: 14

Polynomeilla laskeminen ja kuvaajien piirtäminen

# Polynomilla voi tehdä laskutoimituksia
pol = np.poly1d([1, 4, 3])
print(pol**2)

tulostaa (eksponentit ylärivillä)

    4     3      2
1 x + 8 x + 22 x + 24 x + 9
# Polynomin arvon voi laskea myös useassa pisteessä
X = np.arange(0, 11) # X on array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
Y = pol(X) # Y on array([  3,   8,  15,  24,  35,  48,  63,  80,  99, 120, 143])

# Luodaan polynomin kuvaaja
plt.plot(X, Y)
plt.show()

Kuvaaja on:


Polynomisovitukset

Polynomisovituksia voi tehdä np.polyfit-funktiolla. Tarkastellaan esimerkkiä:

import numpy as np
import matplotlib.pyplot as plt

# Luodaan XY-datat
X = np.array([-4.1, -3.4, -2.2, 1.1, 2.2, 3.3, 4.4, 5.5])
Y = np.array([20.91, 15.46, 8.94, 5.11, 8.94, 14.79, 23.46, 34.15])
# Luodaan raakadatan kuvaaja
plt.plot(X, Y, 'o', label = 'raakadata')

# Sovitusfunktio: np.polyfit(xdata, ydata, polynomin_aste)
# Tehdään raakadatalle toisen asteen polynomisovitus:
kertoimet = np.polyfit(X, Y, 2)
# kertoimet on nyt NumPy-taulukko, joka sisältää sovitetun polynomin kertoimet [a, b, c]:
# array([ 0.9996726 , -0.00626338,  4.00940658])
#         a * x^2      b * x        c * 1

# Tehdään kertoimista sovituspolynomi (np.poly1d)
sovitus = np.poly1d(kertoimet)
# sovitus on nyt poly1d-olio, jolla on sovitetun polynomin kertoimet:
# poly1d([ 0.9996726 , -0.00626338,  4.00940658])
# HUOMAA, että merkintä sovitus[2] palauttaa 2. asteen termin kertoimen 0.9996726!
# np.polyfit- funktion tapauksessa 2. asteen termi on kertoimet[0]!

# Lasketaan sovituspolynomin y-arvot usealle x-arvolle
pol_X = np.linspace(-5.0, 6.0, 50)
pol_Y = sovitus(pol_X)

# Luodaan sovituspolynomin kuvaaja
# Kuvaaja lisätään plt.plot-funktiolla samaan kuvaan raakadatan kanssa
plt.plot(pol_X, pol_Y, color = 'red', label = 'sovitus (2. aste)')
plt.legend(loc = 'upper left', frameon = False)
plt.show()

Lopputulos:


Korrelaatiokerroin

Suorien sovittamisen yhteydessä voi laskea X- ja Y-datojen välisen Pearsonin korrelaatiokertoimen np.corrcoef-funktiolla. Kahden muuttujan välinen korrelaatio kertoo mahdollisesta lineaarisesta riippuvuudesta. Jos korrelaatiokerroin on lähellä arvoa 1, voidaan toisen muuttujan arvo arvioida tietämällä vain toisen arvo. Mitä lähempänä korrelaatiokerroin on nollaa, sitä enemmän arvioituun arvoon liittyy epävarmuutta. 

numpy.corrcoef -funktiolle annetaan X- ja Y-datapisteet, jolloin se palauttaa korrelaatiokertoimet 2×2 taulukkona [[xx, xy], [yx, yy]]. Useimmiten riittää, että poimimme taulukosta xy-korrelaation, eli indeksin [0, 1].

import numpy as np
import matplotlib.pyplot as plt

# Luodaan taulukot mittausarvoista
konsentraatiot = np.array([0.1, 0.2, 0.3, 0.4, 0.5])
absorbanssi = np.array([2.24, 4.02, 6.11, 8.27, 10.56])

# Luodaan kuvaaja raakadatoista
plt.plot(konsentraatiot, absorbanssi, 'o', label = 'raakadata')

# Sovitusfunktio: np.polyfit(xdata, ydata, polynomin_aste)
# Tehdään 1.  asteen polynomisovitus XY-dataan:
kertoimet = np.polyfit(konsentraatiot, absorbanssi, 1)
# kertoimet on nyt NumPy-taulukko, joka sisältää sovitetun polynomin kertoimet

# Tehdään kertoimista polynomi np.poly1d
polynomi = np.poly1d(kertoimet)
# Lasketaan y:n arvot usealle x:n arvolle
pol_X = np.linspace(0.0, 0.6, 60)
pol_Y = polynomi(pol_X)

# Luodaan sovituspolynomin kuvaaja
# Kuvaaja lisätään plot-funktiolla samaan kuvaan raakadatan kanssa
plt.plot(pol_X, pol_Y, color = 'blue', label = 'sovitus (1. aste)')

# Kuvaajan asetukset
plt.xlabel('Konsentraatio (mol/l)')
plt.ylabel('Absorbanssi')
plt.xlim(0, 0.6)
plt.ylim(0, 12)
plt.legend(loc = 'upper left', frameon = False)

plt.show()

# Lasketaan Pearsonin korrelaatiokerroin ja sen neliö
R = np.corrcoef(konsentraatiot, absorbanssi)[0, 1]
R_toiseen = R**2
print(f"Sovituksen R^2 on: {R_toiseen:.3f}")

Lopputuloksena saadaan kuvaaja:


ja tulostus

Sovituksen R^2 on: 0.998

Eli tässä tapauksessa X- ja Y-arvojen välillä on erittäin vahva lineaarinen korrelaatio (absorbanssi on suoraan verrannollinen konsentraatioon).

Tehtävä 4.9.1