SciPy

Kategória: Python.
Alkategória: Python adatszerkezetek.

A SciPy csomag

A SciPy csomag a Scientific Python rövidítése, és - ahogy a nevéből is kitalálható - a természettudományok területét igyekszik lefedni. Ugyanaz a szerzője, mint a NumPy csomagé, számos esetben épít is rá.

Nálam alapból fel volt telepítve; ha mégsem, akkor a szokásos módon telepíthetjük:

pip install scipy

Ez a terület is nagy, igen mély, komoly elmerülést igényel, így itt is tényleg csak a téma felszínét karcoljuk.

Konstansok

A constants modulban igen sok konstans szerepel a természettudományok különbözp területeiről:

from scipy import constants
 
print(constants.Avogadro)      # 6.02214076e+23
print(constants.atmosphere)    # 101325.0
print(constants.pi)            # 3.141592653589793
print(constants.pico)          # 1e-12
print(constants.electron_mass) # 9.1093837015e-31
print(constants.light_year)    # 9460730472580800.0
print(constants.yard)          # 0.9143999999999999
print(constants.acre)          # 4046.8564223999992
print(constants.gibi)          # 1073741824
print(constants.zero_Celsius)  # 273.15
print(constants.horsepower)    # 745.6998715822701

Az összes konstans lekérdezése:

print(dir(constants))
# ['Avogadro', 'Boltzmann', 'Btu', 'Btu_IT', 'Btu_th', 'ConstantWarning',
# 'G', 'Julian_year', 'N_A', 'Planck', 'R', 'Rydberg', 'Stefan_Boltzmann',
# 'Wien', '__all__', '__builtins__', '__cached__', '__doc__', '__file__',
# '__loader__', '__name__', '__package__', '__path__', '__spec__',
# '_obsolete_constants', 'acre', 'alpha', 'angstrom', 'arcmin', 'arcminute',
# 'arcsec', 'arcsecond', 'astronomical_unit', 'atm', 'atmosphere', 'atomic_mass',
# 'atto', 'au', 'bar', 'barrel', 'bbl', 'blob', 'c', 'calorie', 'calorie_IT',
# 'calorie_th', 'carat', 'centi', 'codata', 'constants', 'convert_temperature',
# 'day', 'deci', 'degree', 'degree_Fahrenheit', 'deka', 'dyn', 'dyne', 'e', 'eV',
# 'electron_mass', 'electron_volt', 'elementary_charge', 'epsilon_0', 'erg', 'exa',
# 'exbi', 'femto', 'fermi', 'find', 'fine_structure', 'fluid_ounce', 'fluid_ounce_US',
# 'fluid_ounce_imp', 'foot', 'g', 'gallon', 'gallon_US', 'gallon_imp', 'gas_constant',
# 'gibi', 'giga', 'golden', 'golden_ratio', 'grain', 'gram', 'gravitational_constant',
# 'h', 'hbar', 'hectare', 'hecto', 'horsepower', 'hour', 'hp', 'inch', 'k', 'kgf',
# 'kibi', 'kilo', 'kilogram_force', 'kmh', 'knot', 'lambda2nu', 'lb', 'lbf',
# 'light_year', 'liter', 'litre', 'long_ton', 'm_e', 'm_n', 'm_p', 'm_u', 'mach',
# 'mebi', 'mega', 'metric_ton', 'micro', 'micron', 'mil', 'mile', 'milli', 'minute',
# 'mmHg', 'mph', 'mu_0', 'nano', 'nautical_mile', 'neutron_mass', 'nu2lambda', 'ounce',
# 'oz', 'parsec', 'pebi', 'peta', 'physical_constants', 'pi', 'pico', 'point', 'pound',
# 'pound_force', 'precision', 'proton_mass', 'psi', 'pt', 'short_ton', 'sigma', 'slinch',
# 'slug', 'speed_of_light', 'speed_of_sound', 'stone', 'survey_foot', 'survey_mile',
# 'tebi', 'tera', 'test', 'ton_TNT', 'torr', 'troy_ounce', 'troy_pound', 'u', 'unit',
# 'value', 'week', 'yard', 'year', 'yobi', 'yotta', 'zebi', 'zepto', 'zero_Celsius', 'zetta']

Függvény analízis

A függvények zérushelyeit, minimum és maximum értékeit stb. általános esetben nem egyszerű megtalálni. Számos közelítő módszert dolgoztak ki, melyek közül itt most kettőt vizsgálunk meg.

Függvény gyöke

A lineáris vagy másodfokú függvények nulla helyeit a szokásos módszerekkel, képletekkel ki tudjuk számolni, de mit kezdjünk mondjuk egy $x\cdot\sin x + 4$ függvénnyel? Ezt csak numerikus módszerekkel tudjuk meghatározni, melyhez a SciPy segítséget nyújt:

from scipy.optimize import root
from math import sin
 
my_root = root(lambda x: x * sin(x) + 4, 0)
print(my_root.x)  # [5.46130801]

Itt a root() függvény második paramétere (0) azt jelenti, hogy honnan indulva keresse a minimum helyet. A függvény egyébként számos egyéb paramétert kaphat, többek között 10 módszer közül választhatunk. További részleteket a https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root.html oldalon olvashatunk.

Függvény minimuma

A minimize() függvénnyel kereshetjük meg a függvény minimumpontját ill. az ott felvett értéket:

from scipy.optimize import minimize
 
my_min = minimize(lambda x: 3*x**2 - 2*x + 4, 1)
print(my_min.x)   # [0.33333332]
print(my_min.fun) # 3.666666666666667

Ritka adatok kezelése

Ritka adatok alatt azt értjük, hogy az adatok nagy része nulla. A csr_matrix() segítségével tudunk listát készíteni, ami tartalmazza az x-y koordinátákat és az ott található értéket:

import numpy as np
from scipy.sparse import csr_matrix
 
arr = np.array([
    [0, 0, 0, 0],
    [0, 0, 1, 0],
    [3, 0, 2, 0],
])
print(csr_matrix(arr))
#   (1, 2)    1
#   (2, 0)    3
#   (2, 2)    2

A CSR egyébként a Compressed Sparse Row rövidítése, azaz ez sor orientált eredményt ad. Oszlop orientált eredményt a csc_matrix segítségével kapunk.

Gráf algoritmusok

Egy gráf létrehozása

Hozzuk létre az alábbi gráfot!

Az alábbi példákban az A csúcs 0-val, a B 1-gyel stb. van jelölve.

import numpy as np
from scipy.sparse import csr_matrix
 
adjacency_matrix = np.array([
    [0, 1, 2, 0, 0],
    [1, 0, 0, 5, 0],
    [2, 0, 0, 3, 0],
    [0, 5, 3, 0, 0],
    [0, 0, 0, 0, 0],
])
edge_list = csr_matrix(adjacency_matrix)
print(edge_list)
#   (0, 1)    1
#   (0, 2)    2
#   (1, 0)    1
#   (1, 3)    5
#   (2, 0)    2
#   (2, 3)    3
#   (3, 1)    5
#   (3, 2)    3

Közvetlenül a következő, kissé olvashatatlan módon tudjuk létrehozni:

edge_list = csr_matrix(([1, 2, 1, 5, 2, 3, 5, 3], ([0, 0, 1, 1, 2, 2, 3, 3], [1, 2, 0, 3, 0, 3, 1, 2])), shape=(5, 5))

Összekapcsolt komponensek

Először határozzuk meg az összekapcsolt komponenseket:

from scipy.sparse.csgraph import connected_components
print(connected_components(edge_list))
# (2, array([0, 0, 0, 0, 1]))

Az eredmény értelmezése: számokat rendel a komponensekhez, a példában 0-t és 1-et. Ahogy különálló komponens van, annyi számot hazsnál fel. Az egyes számok azt jelentik, hogy melyik pont melyik komponenshez tartozik. A példában a 0, az 1, a 2 és a 3-as sorszámú pont tartozik ugyanabba, amit 0-val jelölt, a 4-es sorszámú pedig külön, amit 1-gyel jelölt.

Mélységi bejárás

Mélységi bejárás során elindulunk valamerre, és ha már jártunk arra, akkor visszalépünk addig, amíg van valós újabb lehetőség, és arra folytatjuk utunkat. Az eredmény a bejárás sorrendjét adja meg, valamint az, hogy adott csúcshoz honnan érkeztünk:

from scipy.sparse.csgraph import depth_first_order
print(depth_first_order(edge_list, 0))
# (array([0, 1, 3, 2]), array([-9999,     0,     3,     1, -9999]))

A példában a bejárás sorrendje 0, 1, 3, 2 volt. A 0-hoz nem érkeztünk (-9999), az 1-hez a 0-ból érkeztünk, a 2-be a 3-ból érkeztünk, a 3-ba az 1-ből érkeztünk, a 4-be pedig nem értünk oda.

Szélességi bejárás

Szélességi bejárás során először a kiinduló csúcs közvetlen szomszédjait járjuk be, majd azoknak a még be nem jár szomszédait, és ezt folytatjuk mindaddig, amíg tudunk hova menni:

from scipy.sparse.csgraph import breadth_first_order
print(breadth_first_order(edge_list, 0))
# (array([0, 1, 2, 3]), array([-9999,     0,     0,     1, -9999]))

Az eredmény értelmezése: 0, 1, 2, 3 sorrendben jártuk be a csúcsokat. A 0-ba sehonnan sem érkeztünk, mivel onnan indultunk; az 1-be és a 2-be a 0-ból, a 3-ba az 1-ből, a 4-be pedig nem érkeztünk meg. Súlyozatlan gráfban egyébként ezzel a módszerrel találhatjuk meg a legrövidebb utat.

Dijsktra

A Dijkstra algoritmussal olyan irányított súlyozott gráfban keressük meg két csúcs között a legrövidebb utat, melyben nincs negatív súlyú él.

from scipy.sparse.csgraph import dijkstra
print(dijkstra(edge_list, indices=0))
# [ 0.  1.  2.  5. inf]

A SciPy az összes élre kiszámolja a távolságot, és ezzel nem lő túl messzire, mert aszimptotikusan ugyanannyi ideig tart kiszámolni az összeset, mint csak egyet. Az eredmény értelmezése: a 0 sorszámú él 0-ra van a kiindulótól (0), az 1-es sorszámú 1-re, a 2-es 2-re, a 3-as 5-re (0 → 2 → 3, 2+3=5), míg a 4-es elérhetetlen.

Floyd-Warshall

A Floyd-Warshall algoritmus minden csúcspárra kiszámolja a legrövidebb utat:

from scipy.sparse.csgraph import floyd_warshall
print(floyd_warshall(edge_list))
# [[ 0.  1.  2.  5. inf]
#  [ 1.  0.  3.  5. inf]
#  [ 2.  3.  0.  3. inf]
#  [ 5.  5.  3.  0. inf]
#  [inf inf inf inf  0.]]

Az eredmény egy mátrix, amely megadja, hogy melyik csúcsból melyik csúcsba mekkora súlyozott utat kell megtenni.

Bellman-Ford

A Bellman-Ford szintén a legrövidebb utat számolja ki minden csúcspárra, ez viszont negatív élekkel is működik. Negatív összsúlyú kör nem lehet a gráfban.

from scipy.sparse.csgraph import bellman_ford
print(bellman_ford(edge_list))
# [[ 0.  1.  2.  5. inf]
#  [ 1.  0.  3.  5. inf]
#  [ 2.  3.  0.  3. inf]
#  [ 5.  5.  3.  0. inf]
#  [inf inf inf inf  0.]]

Felszíni adatok kezelése

Delaunay-háromszögelés

A módszerrel háromszögekre tudjuk osztani a teret.

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
 
points = np.array([
    [2, 4],
    [3, 4],
    [3, 0],
    [3, 2],
    [2, 2],
    [4, 1]
])
 
triangles = Delaunay(points).simplices
print(triangles)
# [[3 2 5]
#  [3 4 2]
#  [1 3 5]
#  [4 3 0]
#  [3 1 0]]
 
plt.triplot(points[:, 0], points[:, 1], triangles)
plt.scatter(points[:, 0], points[:, 1], color='r')
plt.show()
scipy-delaunay.png

Konvex burok

Konvex burkot a ConvexHull() hívással tudunk létrehozni:

from scipy.spatial import ConvexHull
 
hull = ConvexHull(points)
hull_points = hull.simplices
print(hull_points)
# [[2 5]
#  [1 0]
#  [1 5]
#  [4 0]
#  [4 2]]
plt.scatter(points[:,0], points[:,1])
for simplex in hull_points:
    plt.plot(points[simplex,0], points[simplex,1], 'k-')
plt.show()
scipy-convexhull.png

Távolságok

Sokféle távolság definíció létezik.

Az Euklideszi a közvetlen út hozza:

from scipy.spatial.distance import euclidean
print(euclidean((1, 2), (4, 6))) # 5.0

A Manhatten olyan, hogy egy képzeletbeli négyzetrács oldalain haladhatunk:

from scipy.spatial.distance import cityblock
print(cityblock((1, 2), (4, 6))) # 7

A Hamming azt mutatja meg, hogy a két érték bitjei milyen arányban térnek el:

from scipy.spatial.distance import hamming
print(hamming((True, True, False, False), (True, False, True, False))) # 0.5

Szignifikancia teszt

A szignifikancia teszt számszerűsít eltéréseket.

Egy példán keresztül illusztrálom a lényegét. Tegyük fel, hogy van egy csemegekukorica konzerv, amire az van ráírva, hogy 200 gramm kukorica van benne. Természetesen egyik konzervben sincs atomi pontosságú 200 gramm, attól valamivel el fog térni. Megmérünk grammra pontosan 10 darab konzervet, felírjuk a értékeket; legyenek mondjuk . Tehát 200 körüli értékeket kaptunk. Ebben a példában a null-hipotézis az az, hogy a konzervben 200 gramm kukorica van. A mérések alapján az átlag 197,9 grammra jött ki. A kérdés az, hogy szignifikáns-e ez az eltérés, azaz a gyár csal-e, vagy az eltérés a véletlen műve. Ehhez konkrétan az ún. egy mintás t-tesztet kell használni. A teszt eredménye az ún. p-érték, ami azt mondja meg, hogy mekkora az esélye annak, hogy az eltérés a véletlen műve, feltéve, hogy a null hipotézis igaz. Azaz mekkora annak az esélye, hogy a csemege kukorica mennyisége ténylegesen 197,9 gramm, ha az elméleti értéke 200 gramm. Számoljuk ki!

from scipy.stats import ttest_1samp
 
corn = [195, 198, 202, 201, 192, 197, 199, 196, 199, 200]
print(ttest_1samp(corn, 200).pvalue) # 0.05400623416381284

Az eredmény értelmezése: kb. 5,4% az esélye annak, hogy ez az érték (egészen pontosan: egy legalább ilyen mértékű eltérés) véletlenszerűen is kijöhet. A gyakorlatban előre le szokás rögzíteni egy ún. alfa értéket, melynek tipikus értéke 5%, azaz 0,05, és ha a p-érték az alatti, akkor szignifikánsnak tartjuk az eltérés (jelen esetben kimondjuk, hogy a gyár csal), míg az a feletti értéknél nem tartjuk annak. Ebben a példában eléggé "remeg a léc", de ez alapján még nem mondhatjuk ki a gyár "bűnösségét".

Viszont érdemes megjegyezni, hogy a mintavétel növelésével önmagában lehet csökkenteni a p-értéket, feltéve, ha tényleg "vaj van a gyár füle mögött". Ha még 10 konzervet megvizsgálnánk, és ugyanezek az értékek jönnének ki (azaz az átlag továbbra is 197,9), akkor a p-érték már kb. 0.0045 lesz, azaz fél százalék alatti. Ez esetben már elvethetjük a null-hipotézist, és kimondhatjuk, hogy a konzervben szignifikánsan kevesebb kukorica van, mint amennyit magáról állít.

Számos ilyen teszt létezik, pl. egy eloszlás követ-e egy előre meghatározott mintát (pl. az emberek kb. fele férfi és fele nő; egy előadásra 9 férfi és 6 nő érkezett, szignifikáns-e az eltérés), adott sokaság egy adott eloszlást követ-e stb. Ezeket mind nem nézzük itt meg, mert túlmutat a téma keretein, egyet viszont igen: az ún. több mintás t-tesztet. Ez az egyik leggyakoribb a gyakorlatban, és arra ad választ, hogy az egyik csoportba tartozó elemek szignifikánsan eltérnek-e a másik csoportétól. Tipikus gyakorlati példa: egy gyógyszer hatásának a vizsgálatakor a betegek egyik része ténylegesen megkapja a gyógyszert, míg egy másik része placebót kap. Valahogy számszerűsítjük a javulást (ill. romlást). Nyilván lesznek olyanok is, akik igazi gyógyszert kaptak, és nem javult az állapotuk, és lesznek olyanok is, akiknek javult, mégis placebót kaptak. Tegyük fel, hogy az jön ki eredményül, hogy a tényleges hatóanyagot tartalmazó gyógyszert kapottak eredményei átlagban jobbak, mint a placebót használók esetében. A kérdés az az, hogy az eltérés szignifikáns-e.

A példában normális eloszlású értékeket generálunk:

import numpy as np
from scipy.stats import ttest_ind
 
np.random.seed(1)
print(ttest_ind(np.random.normal(100, 10, 10), np.random.normal(90, 10, 10)).pvalue) # 0.05314
print(ttest_ind(np.random.normal(100, 10, 10), np.random.normal(85, 10, 10)).pvalue) # 0.00017
print(ttest_ind(np.random.normal(100, 8, 10), np.random.normal(90, 8, 10)).pvalue)   # 0.00309
print(ttest_ind(np.random.normal(100, 10, 20), np.random.normal(90, 10, 20)).pvalue) # 0.00024

Folytassuk a fent megkezdett kukoricás példát. Tegyük fel, hogy most a névleges töltőtömeg 100 gramm, és a szórás (a gép bizonytalansága) 10 gramm. Ez kb. azt jelenti, hogy a konzervek kb. kétharmada a 90-110 gramm tartományba fog esni. Tegyük fel, hogy valaki csal a gyárban: úgy értelmezi, hogy ha megengedett a 10 gramm eltérés, akkor azt bele se teszi, hanem konzervenként kilop átlag 10 gramm kukoricát. A cég belső vizsgálatot rendel el, és összehasonlítja két különböző műszak konzerveit. A tesztek és az eredmények az alábbiak:

  • Az első esetben az egyik várható értéke 100, a másiké 90, 10-10 a szórás. 10-10 konzervet vizsgálnak. A p-érték 5% feletti, tehát nem szignifikáns. Bűnügyi hasonlattal élve ez az az eset, amikor tudjuk, hogy a vádlott bűnös (hiszen az egyik várható értéke 100, míg a másiké 90, azaz egyértelműen van különbség), viszont elítélni még nem tudjuk, mert a bizonyítékok még nem elég meggyőzőek.
  • A második esetben növeltük az eltérést 10-ről 15-re. Azaz tegyük fel, hogy nem 10, hanem 15 grammot csalt. Itt már bőven szignifikáns p-értéket eredményezett a teszt. A bűnügyi hasonlattal ez az, hogy a cselekmény súlya olyan, hogy már annyi bizonyítékkal is el lehet ítélni.
  • A harmadik esetben csökkentettük a szórást, azaz a 100-hoz közeli értékek közelebb vannak a 100-hoz, a 90-hez közeliek pedig a 90-hez, mint az első esetben. Azaz picit pontosítottunk a gépek működésén. A csalás itt is egyből szignifikánssá vált.
  • A negyedik esetben pedig ugyanakkora eltérés és szórás mellett több esetet vizsgáltunk. A bűnügyben ez a több bizonyíték gyűjtésének az esete, ami elegendő arra, hogy elítéljék.

Adatok áttekintése

A gyakorlatban sokszor előfordul az, hogy kapunk egy nagy adathalmazt, és ezt se tudjuk, hogy hogyan fogjunk a feldolgozásának. Érdemes először egy áttekintő képet kapni. A {{describe()} függvény ebben segíthet:

from scipy.stats import describe
print(describe(np.random.normal(100, 15, 50)))
# DescribeResult(
#   nobs=50,
#   minmax=(72.13027203329872, 129.506526238821),
#   mean=101.26988392089712,
#   variance=174.4130488987274,
#   skewness=-0.2944156083172776,
#   kurtosis=-0.5333328075803618)

Az alábbi adatokat adja vissza:

  • megfigyelések száma (noobs, number of observations),
  • minimum és maximum értékek,
  • átlag,
  • variancia (azaz a szórás négyzete),
  • torzulás (ez a szimmetria hiányára utal; a negatív érték arra utal, hogy a csúcstól balra több érték van, mint attól jobbra),
  • púposság (a "nyomottságot" jelenti; a negatív érték azt jelenti, hogy a csúcsa lejjebb van a vártnál).

Meglepő egyébként hogy a mediánt (a középső elemet), bizonyos nevezetes centiliseket kihagyták az eredményből. De így is hasznos információkkal szolgáltat.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License