Data-analytiikkaa Pythonilla

Data-analytiikka antaa vastauksia kysymyksiin

Data-analytiikka on tavoitteellista toimintaa: tavoitteena on vastata kysymyksiin. Data-analytiikan avulla vastataan monenlaisiin kysymyksiin:

  • Minkälainen ikäjakauma asiakkaillamme on?
  • Mihin toimintamme osa-alueisiin asiakkaamme ovat tyytymättömiä?
  • Onko asiakkaan iällä yhteyttä asiakastyytyväisyyteen?
  • Miten yrityksen työilmapiiri on muuttunut viimevuodesta?
  • Ketkä asiakkaistamme ovat vaarassa siirtyä kilpailijalle?
  • Keille tuotteen markkinointikampanja kannattaa suunnata?
  • Mikä mainosvaihtoehdoista tehoaa parhaiten kohderyhmään?
  • Mitä oheistuotteita verkkokaupasta ostaneella kannattaa tarjota?
  • Mikä on tuotteen ennustettu kysyntä ensi kuussa?
  • Liittyykö vakuutuskorvaushakemukseen vakuutuspetos?
  • Millä todennäköisyydellä laina-asiakas ei pysty maksamaan lainaansa takaisin?

Data

Tavoitteen (kysymykset, joihin halutaan vastata) asettamisen jälkeen pitää selvittää minkälaista dataa tarvitaan. Data voi olla esimerkiksi:

  • Yrityksen tietokannoista löytyvää dataa (esimerkiksi CRM- ja ERP-järjestelmistä).
  • Erilaisten tiedontuottajien tarjoamaa ilmaista tai maksullista dataa.
  • Varta vasten kyselytutkimuksella tai kokeellisella tutkimuksella kerättyä dataa.
  • Erilaisten sensorien/mittalaitteiden mittaamaa dataa.

Data-analytiikan tarpeisiin datan täytyy olla taulukkomuotoista. Yleisiä data-analytiikkaan sopivia tiedostomuotoja ovat pilkkueroteltu tekstimuoto (.csv) ja Excel-muoto (.xls tai .xlsx). Tietokannoista data haetaan kyselyiden avulla. Nettikyselyohjelmista datan saa yleensä ulos pilkkuerotellussa tekstimuodossa tai Excel-muodossa.

Kun sopiva data on olemassa, niin datasta saadaan vastauksia kysymyksiin seuraavien vaiheiden kautta:

  • Datan valmistelu.
  • Datan kuvailu.
  • Tilastollinen merkitsevyys: Tämä vaihe tulee kyseeseen, jos data pohjautuu isommasta joukosta poimittuun otokseen. Tilastollinen merkitsevyys kertoo, millä varmuudella otoksessa havaittuja eroja ja riippuvuuksia voidaan yleistää isompaan joukkoon.
  • Koneoppiminen ja ennakoiva analytiikka.

Datan valmistelu

Datan valmistelu on yleensä data-analytiikan aikaa vievin vaihe. Ensimmäiseksi kannattaa varmistaa datan taulukkomuotoisuus: muuttujien nimet / kenttien nimet / sarakeotsikot ovat ensimmäisellä rivillä, datassa ei ole tarpeettomia tyhjiä rivejä tai sarakkeita, kuhunkin tilastoyksikköön/havaintoyksikköön liittyvät tiedot ovat yhdellä rivillä. Datan valmistelu voi sisältää muiden muassa seuraavia:

  • Muuttujien uudelleen nimeäminen: jatkotoimet sujuvat sutjakkaammin, jos nimet ovat lyhyitä ja helposti tunnistettavia.
  • Desimaalipilkkujen tarkistaminen: vaikka Suomessa desimaalipilkkuna käytetään pilkkua, niin Pythonissa täytyy käyttää pistettä.
  • Päivämäärien muuntaminen päivämääriksi tunnistettavaan muotoon.
  • Mittayksiköiden tarkistaminen ja tarvittavien muunnosten tekeminen.
  • Puuttuvien arvojen käsittely: poistetaanko puuttuvia arvoja sisältävät rivit, korvataanko puuttuvat arvot jollain, miten puuttuvia arvoja merkitään?
  • Uusien muuttujien laskeminen: esimerkiksi summamuuttuja useasta mielipidemuuttujasta, tilauksen hinta tilausmäärän ja yksikköhinnan avulla jne.
  • Arvojen luokittelu ja uudelleenkoodaaminen: esimerkiksi ikäluokat iän arvoista.

Datan kuvailu

Datan kuvailu voi sisältää seuraavia:

  • Lukumäärä- ja prosenttiyhteenvetojen laskeminen (frekvenssitaulukot).
  • Tilastollisten tunnuslukujen laskeminen (keskiarvo, keskihajonta, viiden luvun yhteenveto).
  • Riippuvuuksien tarkastelu (ristiintaulukoinnit, korrelaatiot).
  • Prosenttimuutosten laskeminen aikasarjoille.
  • Aikasarjojen tarkastelu viivakuvioina.
  • Liukuvien keskiarvojen esittäminen aikasarjojen yhteydessä.

Kuvailun tuloksia kannattaa visualisoida ja havainnollistaa hyvin viimeistellyillä taulukoilla ja kuvioilla.

Tilastollinen merkitsevyys

Jos käytetty data on otos isommasta perusjoukosta, niin tulokset kuvaavat otosta. Jos tarkoituksena on arvioida koko perusjoukkoa, niin otoksessa havaittujen erojen ja riippuvuuksien tilastollinen merkitsevyys kertoo, millä varmuudella eroja ja riippuvuuksia voidaan yleistää perusjoukkoon.

Koneoppiminen ja ennakoiva analytiikka

Koneoppimisen malleilla voidaan luokitella (asiakkaat luottoriski-asiakkaisiin ja muihin, vakuutuskorvaushakemukset selviin tapauksiin ja petokselta haiskahtaviin jne.) ja ennakoida (tulevaa kysyntää jne.). Koneoppiminen perustuu siihen, että kone oppii käytettävän mallin parametrit olemassa olevasta datasta ja tämän jälkeen mallia voidaan soveltaa uuteen dataan.

Koneoppimisalgoritmit voidaan luokitella  seuraavasti (suomennokset eivät ole vakiintuneita):

  • Supervised learning (ohjattu oppiminen): Algoritmi opetetaan opetusdatalla (training data). Esimerkiksi roskapostisuodatin opetetaan sähköpostidatalla, jossa on erilaisia tietoja kustakin sähköpostiviestistä sekä tieto siitä oliko sähköpostiviesti roskapostia. Tämän datan perusteella muodostuu malli, jota käyttäen tulevista sähköpostiviesteistä voidaan tunnistaa roskapostiviestit.
  • Unsupervised learning (ohjaamaton oppiminen): Esimerkiksi asiakkaiden jakaminen asiakassegmentteihin.
  • Reinforcement learning (vahvistusoppiminen): Algoritmi suorittaa toimia ja saa niistä palautetta palkkioiden ja rangaistuksen muodoissa. Algoritmi oppii saamistaan palkkioista ja rangaistuksista. Vahvistettua oppimista käytetään esimerkiksi robotiikassa.

Seuraavassa jaotellaan ohjattu ja ohjaamaton oppiminen edelleen alatyyppeihin:

kone1

Ohjattu oppiminen

Label tarkoittaa ennakoitavaa asiaa (selitettävää muuttujaa).

Diskreetti label

Jos ennakoitava asia on diskreetti (epäjatkuva), niin kyseeseen tulevat luokittelua suorittavat algoritmit, esimerkiksi logistinen regressio tai päätöspuut.

Esimerkkejä, joissa on diskreetti label:

  • Roskapostisuodatin: Label on tieto siitä, onko sähköpostiviesti roskapostia vai ei?
  • Lääketieteellinen diagnoosi: Label on tieto siitä, onko tutkitulla potilaalla tietty sairaus vai ei?
  • Vakuutuspetosten tunnistaminen: Label on tieto siitä, liittyykö korvaushakemukseen petos vai ei?

Jatkuva label

Jos ennakoitava asia on jatkuva, niin kyseeseen tulevat esimerkiksi regressiomallit ja aikasarjaennustamisen menetelmät. Esimerkkejä, joissa on jatkuva label:

  • Osakehuoneiston hinnan ennustaminen: Label on asunnon hinta.
  • Kysynnän ennustaminen aikaisemman kysynnän perusteella: Label on kysyntä.

Ohjaamaton oppiminen

Ohjaamattomassa oppimisessa ei ole labelia (selitettävää/ennakoitavaa muuttujaa). Ohjaamattoman oppimisen algoritmi muodostaa mallin suoraan datasta. Esimerkkinä asiakassegmenttien määrittäminen asiakasdatan pohjalta. Paljon käytetty algoritmi on k-means clustering.

Jos datassa on paljon muuttujia, jotka mittaavat osittain samoja asioita, niin dimensionality reduction -tyyppisillä algoritmeilla voidaan pienentää muuttujien määrää yhdistämällä niitä kimpuiksi. Tunnetuin algoritmi tähän tarkoitukseen on pääkomponenttianalyysi.

Data-analytiikkaa Pythonilla

Jos aiot käyttää Pythonia data-analytiikassa, niin kannattaa aloittaa asentamalla Anaconda.

Mainokset

Kuviot Pythonilla 3

Kategoristen muuttujien jakauman esitetän lukumäärä- tai prosenttipylväitä käyttäen. Näistä voit lukea lisää artikkeleistani Kuviot Pythonilla 1 ja Kuviot Pythonilla 2. Tässä artikkelissa esitän esimerkkejä määrällisen muuttujien kuvaamisesta. Tämän artikkelin ohjelmakoodin ja tulosteet löydät GitHubista:

https://github.com/taanila/tilastoapu/blob/master/kuviot3.ipynb

Jos kopioit koodia itsellesi, niin kannattaa käyttää GitHubia. Tästä artikkelista kopioidut koodit eivät välttämättä toimi oikein.

Oletan, että lukijalla on asennettuna Anaconda ja sen mukana tuleva Jupyter notebook.

Ensimmäiseksi minun täytyy ottaa käyttöön tämän artikkelin kuvioissa tarvittavat kirjastot:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
plt.style.use('seaborn-white')

Tyylin ’seaborn-white’ sijasta voit käyttää sinulle mieluista tyyliä.

Avaan seaborn-kirjastoon sisältyvän esimerkkiaineiston tips (tietoja ravintola-asiakkaiden ravintolakäynneistä):

tips = sns.load_dataset('tips')
tips.head()

tips

Histogrammi

Ensimmäiseksi kuvaan ravintola-asiakkaiden laskujen loppusumman jakauman histogrammina:

ax=tips['total_bill'].plot.hist(bins=5)
ax.set_xlabel('total_bill')

hist1

Lisämääre bins=5 määrittää luokkien (histogrammin pylväiden) lukumäärän.

Voin myös itse määrittää luokkarajat:

ax=tips['total_bill'].plot.hist(bins=[0,10,20,30,40,50,60], 
   edgecolor='white')
ax.set_xlabel('total_bill')

hist2

Antamani luokkarajat tuottavat luokat [0,10), [10,20), [20,30), [30,40), [40,50) ja [50,60]. Luokan alaraja siis sisältyy aina luokkaan, mutta luokan yläraja ei sisälly luokkaan viimeistä luokkaa lukuun ottamatta.

Lisämääre edgecolor=’white’ lisää pylväille valkoisen reunan.

Seuraavaksi kuvaan palvelurahan määrän jakaumaa. Jos haluan esittää histogrammissa lukumäärien sijasta prosentit, niin voin käyttää esimerkiksi seuraavaa esoteerisen näköistä koodia:

weights=np.ones_like(tips['tip'])/float(len(tips['tip']))
ax=tips['tip'].plot.hist(bins=[0,2,4,6,8,10], weights=weights, 
   edgecolor='white')
ax.set_xlabel('tip')
ax.set_ylabel('%')
vals = ax.get_yticks()
ax.set_yticklabels(['{:.0f} %'.format(y*100) for y in vals])

hist3

weights vaatinee selitystä:

Otin edellä käyttöön numpy-kirjaston (import numpy as np). Numpy (numerical python) on peruskirjasto numeerisen tiedon käsittelyyn. Tässä käytän np.ones_like toimintoa, joka korvaa kunkin ravintola-asiakkaan palvelurahan ykkösellä. Jakamalla ykkösen palvelurahaa maksaneiden lukumäärällä (float(len(tips[’tip’]) saan selville prosenttiosuuden kaikista palvelurahaa maksaneista. Kaikkiaan weights=np.ones_like(tips[’tip’])/float(len(tips[’tip’])) muodostaa sarjan prosenttiosuuksia, joilla voin painottaa histogrammin lukumääriä ja näin lopulta saan histogrammiin prosentit.

Todennäköisyysjakauman estimointi

Voin edetä histogrammista taustalla olevan todennäköisyysjakauman tiheysfunktion estimointiin käyttämällä seaborn-kirjaston distplot-toimintoa:

sns.distplot(a=tips['tip'].dropna(), bins=[0,2,4,6,8,10])

distplot1

distplot antaa virheilmoituksen puuttuvista arvoista. Tämän vuoksi poistan dropna()-funktiolla mahdollisia puuttuvia arvoja sisältävät rivit. Kuvion pystyakselin arvot ovat todennäköisyysjakauman tiheysfunktion arvoja. Todennäköisyysjakauma estimoidaan ei-parametrisella KDE-menetelmällä (kernel density estimation). Voit lukea Wikipediasta lisää KDE-menetelmästä https://en.wikipedia.org/wiki/Kernel_density_estimation

Lisämääreellä kde=False voin laatia pelkän histogrammin:

sns.distplot(a=tips['tip'].dropna(), bins=[0,2,4,6,8,10], 
   kde=False, rug=True)

distplot2

Luokkien lukumäärien lisäksi voin tarkastella arvojen jakaumia luokkien sisällä alareunan viivoista. Nämä tulivat kuvioon lisämääreen rug=True ansiosta.

Ruutu- ja janakaavio

Ruutu- ja janakaavio on tehokas tapa määrällisen muuttujan jakauman tarkasteluun kategorisen muuttujan määräämissä luokissa. Voin esimerkiksi tarkstella palvelurahan suuruutta eri kokoisten seurueiden kohdalla:

ax=tips.boxplot(column='tip', by='size')
ax.set_title('')
plt.suptitle('')
ax.set_ylabel('tip')

box1

Ruutu- ja janakaavio vaatii lähtötietoinaan kuvattavan muuttujan (column) ja luokittelevan muuttujan (by). Jos ruutu- ja janakaavio ei ole sinulle entuudestaan tuttu, niin lue lisää artikkelistani Ruutu- ja janakaavio.

Voin laatia ruutu- ja janakaavion kätevämmin seaborn-kirjaston avulla. Esimerkiksi palvelurahat seurueen koon mukaan luokiteltuna sukupuolittain:

sns.boxplot(x='size', y='tip', hue='sex', data=tips)

box2

Pairplot

Kahden muuttujan välisiä riippuvuuksia voin tarkatella seaborn-kirjaston pairplot-kuviona.

Käytän esimerkkinä aineistoa, jossa on eri aiheiden osaamista/lahjakkuutta kuvaavia pistemääriä ja arvosanoja sekä opintomenestystä kuvaava pistemäärä.

opintomenestys = pd.read_excel('http://taanila.fi/
   opintomenestys.xlsx')
opintomenestys.head()

opintomenestys
Voin laatia kaikki parittaiset hajontakuviot useammasta muuttujasta:

sns.pairplot(opintomenestys[['verbaalinen','looginen','kielet',
   'matematiikka','opintomenestys']], kind='reg')

pairplot1

Yllä on näkyvillä vain osa kuvioista. Tuloksena saan jokaisesta muuttujasta histogrammin ja kaikki parittaiset hajontakaaviot pienimmän neliösumman regressiosuoralla täydennettynä. Voin tarkastella asiaa myös sukupuolittain:

sns.pairplot(opintomenestys[['verbaalinen','looginen','kielet',
   'matematiikka','opintomenestys','sukupuoli']], 
   hue='sukupuoli', kind='reg')

pairplot2

Histogrammien sijasta saan KDE-menetelmällä estimoidut todennäköisyysjakaumat miehille (oranssi) ja naisille (sininen). Jos haluan histogrammit, niin voin käyttää lisäparametria diag_kind=’hist’.

Jointplot

Yksittäisen muuttujaparin hajontakuvion saan kätevimmin seaborn-kirjaston jointplot-toiminnolla:

sns.jointplot(x='kielet', y='opintomenestys', data=opintomenestys, 
   kind='reg')

jointplot1

Kuvion reunoilla on muuttujien histogrammit sekä KDE-menetelmällä estimoidut todennäköisyysjakumat.

Tarvittaessa voin laskea kuvioon korrelaatiokertoimen ja sen merkitsevyyttä mittaavan p-arvon. Korrelaatiokertoimen laskemiseen tarvittavat funktiot löytyvät scipy.stats-kirjastosta. Pearsonin korrelaatiokerroin löytyy pearsonr-nimellä ja Spearmanin järjestyskorrelaatiokerroin spearmanr-nimellä.

from scipy.stats import spearmanr
sns.jointplot(x='matematiikka', y='opintomenestys', data=opintomenestys, 
   kind='reg', stat_func=spearmanr)

jointplot2

Aikasarjaennustaminen 3

Tämä artikkeli on jatkoa artikkeleille Aikasarjaennustaminen 1 ja Aikasarjaennustaminen 2.

Edellisen artikkelin Aikasarjaennustaminen 2 lopussa totesin, että esimerkkinä käyttämässäni aikasarjassa on neljän vuosineljänneksen välein toistuvaa kausivaihtelua, joka on syytä huomioida ennustamisessa. Tässä artikkelissa tarkastelen kausivaihtelun huomioivaa Holt-Winterin menetelmää.

Holt-Winterin tulomallissa aikasarjan tason L (level) hetkellä t määrittää lauseke

Lt = alfa * Yt/St-s + (1 – alfa)(Lt-1 + Tt-1)

Yllä Yt on viimeisin havainto, St-s on edellisen vastaavan periodin kausivaihtelu ja  Tt-1 on edellinen trendi.

Trendille T hetkellä t saadaan arvio lausekkeesta

Tt = beta * (Lt – Lt-1) + (1 – beta) * Tt-1

Kausivaihtelulle S hetkellä t saadaan arvio lausekkeesta

St = gamma * Yt/Lt + (1 – gamma) * St-s

Ennuste hetkelle t + p saadaan

(Lt + pTt)St-s

Yllä on kyse Holt-Winterin tulomallista, jossa kausivaihtelu huomioidaan kausivaihtelukertoimena. Holt-Winterin mallia voidaan soveltaa myös summamallina, jolloin kausivaihtelu huomioidaan lisättävänä kausivaihteluterminä. Tulomalli soveltuu paremmin tilanteisiin, joissa kausivaihtelukomponentin suuruus vaihtelee aikasarjan tason L mukaan. Summamalli soveltuu tilanteisiin, joissa kausivaihtelukomponentin suuruus ei riipu aikasarjan tasosta L.

Mallin parametrit alfa, beta ja gamma pyritään määrittämään siten että virheiden keskihajonta RMSE saa pienimmän mahdollisen arvonsa.

Python toteutus

Tämän artikkelin ohjelmakoodin ja tulosteet löydät GitHubista:

https://github.com/taanila/tilastoapu/blob/master/forecast3.ipynb

Jos kopioit koodia itsellesi, niin kannattaa käyttää GitHubia. Tästä artikkelista kopioidut koodit eivät välttämättä toimi oikein (esimerkiksi sisennykset voivat kopioitua virheellisesti).

Oletan, että lukijalla on asennettuna Anaconda ja sen mukana tuleva Jupyter notebook.

Kirjastot

Tarvitsen monia ohjelmakirjastoja: pandas dataframen käsittelyyn, matplotlib kuvioiden tekemiseen, statsmodels eksponentiaaliseen tasoitukseen, sklearn ennustevirhettä kuvaavien tunnuslukujen laskemiseen ja math neliöjuuren laskemiseen.

import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
from statsmodels.tsa.api import Holt
from sklearn.metrics import mean_squared_error, mean_absolute_error
from math import sqrt

Aikasarja

Tallennan vuosineljänneksittäisiä myyntilukuja sisältävän aikasarjan data-nimiseen dataframeen:

series=[500,350,250,400,450,350,200,300,350,200,150,
   400,550,350,250,550,550, 400,350,600,750,500,400,650,850]
index=pd.date_range('2000-03-31', periods=25, freq='Q')
data=pd.DataFrame(series, index=index)
data.columns=['Demand']

Ennustemalli

Tämän jälkeen sovitan eksponentiaalisen mallin aikasarjaan ja tallennan mallin mukaiset sovitetut arvot datan Predict-sarakkeeseen. Jos haluan käyttää summamallia, niin korvaan parametrin seasonal=’mul’ parametrilla seasonal=’add’.

fit1 = ExponentialSmoothing(data['Demand'], seasonal_periods=4, 
   trend='add', seasonal='mul').fit()
data['Predict']=fit1.fittedvalues 
data.plot()

holtwinter1

Silmämäärin katsottuna Holt-Winterin malli istuu hyvin aikasarjaan. Tarkistan asian virheiden suuruutta kuvaavista tunnusluvuista. Ensiksi RMSE:

sqrt(mean_squared_error(data['Demand'], data['Predict']))

Tulokseksi saan RMSE = 53 ( yksinkertaisessa eksponentiaalisessa tasoituksessa 150)

MAD:

mean_absolute_error(data['Demand'], data['Predict'])

Tulokseksi saan MAD = 42 (yksinkertaisessa eksponentiaalisessa tasoituksessa 128)

Virheitä kannattaa tarkastella myös sellaisenaan:

data['Resid']=fit1.resid
data['Resid'].plot()

holtwinter2

Virheiden vaihtelu näyttää satunnaiselta niin kuin hyvässä mallissa pitää ollakin.

Ennuste

Lasken vielä ennusteet neljälle seuraavalle neljännekselle omaan dataframeen ja piirrän samaan kuvioon alkuperäisen aikasarjan ja neljän  neljänneksen ennusteet.

index=pd.date_range('2006-06-30', periods=4, freq='Q')
datap=pd.DataFrame(fit1.forecast(4), index=index)
datap.columns=['Predict']
data['Demand'].plot()
datap['Predict'].plot()

holtwinter3

Ennusteestakin (oranssi) näkee, että kausivaihtelu on mallissa huomioitu.

Mallin parametreista selviävät alfan (smoothing_level), betan (smoothing_slope), ja gamman (smoothing_seasonal) arvot sekä tason, trendin ja kausivaihtelun alkuarvot :

fit1.params

{’damping_slope’: nan,
’initial_level’: 571.4286398129714,
’initial_seasons’: array([0.84785489, 0.56435295, 0.41577146, 0.71556852]),
’initial_slope’: 18.700003064693195,
’lamda’: None,
’remove_bias’: False,
’smoothing_level’: 0.8700136796734163,
’smoothing_seasonal’: 0.0,
’smoothing_slope’: 1.818780385451967e-19,
’use_boxcox’: False}

Aikasarjaennustaminen 2

Tämä artikkeli on jatkoa yksinkertaista eksponentiaalista tasoitusta käsittelevälle artikkelille Aikasarjaennustaminen 1. Tässä artikkelissa käytän kaksinkertaista eksponentiaalista tasoitusta eli Holtin mallia, joka huomioi myös trendin.

Holtin mallissa aikasarjan tason L (level) hetkellä t määrittää lauseke

Lt = alfa * Yt + (1 – alfa) * (Lt-1 + Tt-1)

Yllä Yt on viimeisin havainto ja Tt-1 on edellinen trendi. Trendille hetkellä t saadaan arvio lausekkeesta

Tt = beta * (Lt – Lt-1) + (1 – beta) * Tt-1

Ennuste hetkelle t+p saadaan

Lt + pTt

Mallin parametrit alfa ja beta pyritään määrittämään siten että virheiden keskihajonta RMSE saa pienimmän mahdollisen arvonsa.

Python toteutus

Tämän artikkelin ohjelmakoodin ja tulosteet löydät GitHubista:

https://github.com/taanila/tilastoapu/blob/master/forecast2.ipynb

Jos kopioit koodia itsellesi, niin kannattaa käyttää GitHubia. Tästä artikkelista kopioidut koodit eivät välttämättä toimi oikein (esimerkiksi sisennykset voivat kopioitua virheellisesti).

Oletan, että lukijalla on asennettuna Anaconda ja sen mukana tuleva Jupyter notebook.

Kirjastot

Tarvitsen monia ohjelmakirjastoja: pandas dataframen käsittelyyn, matplotlib kuvioiden tekemiseen, statsmodels eksponentiaaliseen tasoitukseen, sklearn ennustevirhettä kuvaavien tunnuslukujen laskemiseen ja math neliöjuuren laskemiseen.

import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
from statsmodels.tsa.api import Holt
from sklearn.metrics import mean_squared_error, mean_absolute_error
from math import sqrt

Aikasarja

Tallennan vuosineljänneksittäisiä myyntilukuja sisältävän aikasarjan data-nimiseen dataframeen:

series=[500,350,250,400,450,350,200,300,350,200,150,
   400,550,350,250,550,550,400,350,600,750,500,400,650,850]
index=pd.date_range('2000-03-31', periods=25, freq='Q')
data=pd.DataFrame(series, index=index)
data.columns=['Demand']

Ennustemalli

Tämän jälkeen sovitan eksponentiaalisen mallin aikasarjaan ja tallennan mallin mukaiset sovitetut arvot datan Predict-sarakkeeseen. Viivakuviosta voin arvioida silmämäärin mallin sopivuutta havaintoihin.

fit1 = Holt(data['Demand']).fit()
data['Predict']=fit1.fittedvalues
data.plot()

holt1

Silmämäärin katsottuna Holtin malli ei vaikuta juurikaan paremmalta kuin yksinkertainen eksponentiaalinen tasoitus. Tarkistan asian vielä virheiden suuruutta kuvaavista tunnusluvuista. Ensiksi RMSE:

sqrt(mean_squared_error(data['Demand'], data['Predict']))

Tulokseksi saan RMSE = 147 ( yksinkertaisessa eksponentiaalisessa tasoituksessa 150)

MAD:

mean_absolute_error(data['Demand'], data['Predict'])

Tulokseksi saan MAD = 125 (yksinkertaisessa eksponentiaalisessa tasoituksessa 128)

Tunnuslukujenkaan valossa malli ei vaikuta paljoa paremmalta.

Ennuste

Lasken vielä ennusteet neljälle seuraavalle neljännekselle omaan dataframeen ja piirrän samaan kuvioon alkuperäisen aikasarjan ja neljän neljänneksen ennusteet.

index=pd.date_range('2006-06-30', periods=4, freq='Q')
datap=pd.DataFrame(fit1.forecast(4), index=index)
datap.columns=['Predict']
data['Demand'].plot()
datap['Predict'].plot()

holt3

Se, että ennusteessa on huomioitu trendi, näkyy kuviossa hyvin.

Mallin parametreista selviävät alfan (smoothing_level) ja betan (smoothing_slope) arvot:

fit1.params

{’damping_slope’: nan,
’initial_level’: 409.18761409566923,
’initial_seasons’: array([], dtype=float64),
’initial_slope’: 0.0,
’lamda’: None,
’remove_bias’: False,
’smoothing_level’: 0.25213126982732664,
’smoothing_seasonal’: nan,
’smoothing_slope’: 0.2521311899925906,
’use_boxcox’: False}

Aikasarjan tasoon liittyvä alfa on noin 0,25 ja trendiin liittyvä beta on noin 0,25.

Aikasarjan lähempää tarkastelua

Löytääkseni paremman ennustemallin tarkastelen aikasarjaa hieman lähemmin purkamalla sen komponentteihin.

from statsmodels.tsa.api import seasonal_decompose
seasonal_decompose(data['Demand']).plot()

Tuloksena saan neljä kuviota:

  • alkuperäinen aikasarja
  • aikasarjasta erotettu trendi
  • aikasarjan kausivaihtelu
  • aikasarjan jäljelle jäänyt osa trendin ja kausivaihtelun poistamisen jälkeen.

holt4

Aikasarjassa on erotettavissa selkeä neljän vuosineljänneksen jaksoissa toistuva kausivaihtelu, jota ennustemallini ei huomioinut. Asiaa voin tarkastella myös autokorrelaatioiden avulla. Autokorrelaatio tarkoittaa aikasarjan korrelaatiota viivästetyn aikasarjan kanssa, esimerkiksi aikasarjan korrelaatio neljän vuosineljänneksen takaisiin aikasarjan arvoihin. Autokorrelaatio voidaan laskea eri viiveille. Tämän voin tehdä pandas-kirjaston autocorrelation_plot-toiminnolla:

from pandas.plotting import autocorrelation_plot
autocorrelation_plot(data['Demand'])

holt5

Vaaka-akselilla on viive (lag) ja pystyakselilla autokorrelaatiokertoimen arvo. Huomaan, että viiveen 4 kohdalla on suurehko korrelaatio. Tämä viittaa 4 vuosineljänneksen mittaiseen kausivaihtelujaksoon. Kuvion katkoviivat edustavat tilastollisesti merkitsevän korrelaation rajoja. Viiveen 4 kohdalla korrelaatio on katkoviivan yläpuolella ja näin ollen tilastollisesti merkitsevä.

Seuraavassa artikkelissa Aikasarjaennustaminen 3 laadin ennustemallin, jossa myös 4 vuosineljänneksen jaksoissa toistuva kausivaihtelu on huomioitu.

 

Aikasarjaennustaminen 1

Aikasarjaennustamisessa oletan, että toteutuneiden havaintojen muodostama aikasarja sisältää informaatiota, joka auttaa tulevien havaintojen ennustamisessa. Ennustusmenetelmä riippuu siitä, minkälaista systemaattista vaihtelua aikasarjassa esiintyy. Eksponentiaalisia tasoitusmenetelmiä käytettäessä on kolme päävaihtoehtoa:

  • Yksinkertainen eksponentiaalinen tasoitus aikasarjoille, joissa ei ole trendiä eikä kausivaihtelua.
  • Kaksinkertainen eksponentiaalinen tasoitus eli Holtin menetelmä aikasarjoille, joissa on trendi, mutta ei kausivaihtelua.
  • Kolminkertainen eksponentiaalinen tasoitus eli Holt-Winterin menetelmä aikasarjoille, joissa on sekä trendi että kausivaihtelu.

Tämä artikkeli käsittelee yksinkertaista eksponentiaalista tasoitusta. Yksinkertaisessa eksponentiaalisessa tasoituksessa ennuste lasketaan seuraavasti:

alfa*edellinen havainto + (1 – alfa)*edellinen ennuste

Ennuste saadaan viimeisimmän havainnon ja siihen liittyneen ennusteen painotettuna summana. Painokerroin alfa on välillä 0 – 1 oleva luku, joka ilmaisee, kuinka suurella painolla edellistä havaintoa painotetaan ennustetta laskettaessa:

    • Jos alfa on 0, niin ennuste on sama kuin edellinen ennuste.
    • Jos alfa on 1, niin ennuste on sama kuin edellinen havainto.
    • Suuret alfan arvot antavat ennusteita, jotka reagoivat herkästi aikasarjassa esiintyvään vaihteluun, koska viimeisimmillä havainnoilla on suurempi paino.
    • Pienet alfan arvot tasoittavat voimakkaasti aikasarjan vaihtelua.

 

Alfan arvo valitaan yleensä siten että keskimääräinen ennustevirhe saadaan mahdollisimman pieneksi. Voin kirjoittaa ennusteen laskentakaavan myös muotoon

edellinen ennuste + alfa*(edellinen havainto – edellinen ennuste)

Ennustetta siis korjataan jokaisen toteutuneen havainnon jälkeen korjaustermillä alfa*edellisen ennusteen virhe.

Python toteutus

Tämän artikkelin ohjelmakoodin ja tulosteet löydät GitHubista:

https://github.com/taanila/tilastoapu/blob/master/forecast1.ipynb

Jos kopioit koodia itsellesi, niin kannattaa käyttää GitHubia. Tästä artikkelista kopioidut koodit eivät välttämättä toimi oikein (esimerkiksi sisennykset voivat kopioitua virheellisesti).

Oletan, että lukijalla on asennettuna Anaconda ja sen mukana tuleva Jupyter notebook.

Kirjastot

Tarvitsen monia ohjelmakirjastoja: pandas dataframen käsittelyyn, matplotlib kuvioiden tekemiseen, statsmodels eksponentiaaliseen tasoitukseen, sklearn ennustevirhettä kuvaavien tunnuslukujen laskemiseen ja math neliöjuuren laskemiseen.

import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
from statsmodels.tsa.api import SimpleExpSmoothing
from sklearn.metrics import mean_squared_error, mean_absolute_error
from math import sqrt

Aikasarja

Tallennan vuosineljänneksittäisiä myyntilukuja sisältävän aikasarjan data-nimiseen dataframeen:

series=[500,350,250,400,450,350,200,300,350,200,150,400,
   550,350,250,550,550,400,350,600,750,500,400,650,850]
index=pd.date_range('2000-03-31', periods=25, freq='Q')
data=pd.DataFrame(series, index=index)
data.columns=['Demand']

Huomaa, kuinka kätevästi määrittelen vuosineljännekset pandas-kirjaston date_range-toiminnolla.

Ennustemalli

Sovitan eksponentiaalisen mallin aikasarjaan ja tallennan mallin mukaiset sovitetut arvot datan Predict-sarakkeeseen. Viivakuviosta voin arvioida silmämäärin mallin sopivuutta havaintoihin.

fit1 = SimpleExpSmoothing(data['Demand']).fit()
data['Predict']=fit1.fittedvalues
data.plot()

exp1

Lukiessasi seuraavia artikkeleitani, huomaat, että parempiakin ennustemalleja on tälle aikasarjalle löydettävissä.

Mallin sopivuutta voin mitata laskemalla ennustevirheiden (ennuste-toteutunut) neliöiden keskiarvon ja ottamalla siitä neliöjuuren (RMSE=root mean square error). Tämä kuvaa ennustevirheiden keskihajontaa.

sqrt(mean_squared_error(data['Demand'], data['Predict']))

Tulokseksi saan noin 150. Tätä kannattaa verrata seuraavissa artikkeleissa laskettuihin vastaaviin lukuihin Holtin ja Holt-Winterin malleissa. Karkeasti voisi todeta: Mitä pienempi RMSE sitä parempi malli.

Toinen usein käytetty mallin sopivuutta mittaava luku on virheiden itseisarvojen keskiarvo (MAD=mean absolute deviation):

mean_absolute_error(data['Demand'], data['Predict'])

Tämän arvoksi saan noin 128.

Virheitä kannattaa tarkastella myös sellaisenaan:

data['Resid']=fit1.resid
data['Resid'].plot()

exp2

Hyvässä ennustemallissa virheissä ei saisi olla trendiä tai systemaattista vaihtelua. Tässä virheen vaihtelu näyttää melko systemaattiselta, joten ennustemalli ei ole paras mahdollinen.

Ennuste

Lasken vielä ennusteet neljälle seuraavalle neljännekselle omaan dataframeen ja piirrän samaan kuvioon alkuperäisen aikasarjan ja neljän neljänneksen ennusteet. Yksinkertainen eksponentiaalinen tasoitus antaa tuleville ajankohdille vakioennusteen, jota toki voidaan korjata jokaisen uuden havainnon jälkeen.

index=pd.date_range('2006-06-30', periods=4, freq='Q')
datap=pd.DataFrame(fit1.forecast(4), index=index)
datap.columns=['Predict']
data['Demand'].plot()
datap['Predict'].plot()

exp3

Mallin parametreihin voin tutustua tarkemmin komennolla

fit1.params

{’damping_slope’: nan,
’initial_level’: 388.1086015386066,
’initial_seasons’: array([], dtype=float64),
’initial_slope’: nan,
’lamda’: None,
’remove_bias’: False,
’smoothing_level’: 0.32274949027017524,
’smoothing_seasonal’: nan,
’smoothing_slope’: nan,
’use_boxcox’: False}

Tästä selviää, että mallissa käytetty alfa on noin 0,32. Tämä alfa on määräytynyt siten että RMSE saa pienimmän mahdollisen arvon. Mallin sovitusvaiheessa voin määrittää alfan myös itse, esimerkiksi:

fit1 = SimpleExpSmoothing(data['Demand']).fit
   (smoothing_level=0.6)

Seuraavaksi

Esimerkkinä käyttämässäni aikasarjassa on melko helppo erottaa alkupään laskeva trendi ja loppupään nouseva trendi. Seuraavassa artikkelissani Aikasarjaennustaminen 2 yritän sovittaa aikasarjaan mallin, joka huomioi trendin.

Kuviot Pythonilla 1

Tässä artikkelissa perehdyn kuvioiden laatimisen (pylväskuviot) alkeisiin. Tämän artikkelin ohjelmakoodin ja tulosteet löydät GitHubista:

https://github.com/taanila/tilastoapu/blob/master/kuviot1.ipynb

Jos kopioit koodia itsellesi, niin kannattaa käyttää GitHubia. Tästä artikkelista kopioidut koodit eivät välttämättä toimi oikein.

Oletan, että lukijalla on asennettuna Anaconda ja sen mukana tuleva Jupyter notebook.

Pandas-kirjaston dataframe-tietorakenne on keskeinen datojen analysoinnissa: avaan analysoitavan datan dataframeen ja lasketut tulostaulukot ovat enimmäkseen dataframe-muotoisia. Tässä ja muissa kuvioita käsittelevissä artikkeleissani tarkastelen kuvioiden luomista dataframe-muotoisesta datasta.

Matplotlib

Kuvioiden luomisen kulmakivi on matplotlib-kirjasto. Yleensä riittää, että otan käyttöön matplotlib.pyplot-kirjaston (ja tietenkin pandas-kirjaston):

import pandas as pd
import matplotlib.pyplot as plt

Jos käytän Jupyter-notebookia, niin suoritan myös komennon:

%matplotlib inline

Tämän ansiosta kuviot tulostuvat Jupyter-notebookiin ilman erillistä tulostuskomentoa.

Yksinkertainen pylväskuvio

Ensimmäisenä esimerkkinä kuvaan yksinkertaisen myyntilukuja sisältävän taulukon pylväskuviona. Luon myynnit-nimisen dataframen, jossa on neljän eri alueen myyntiluvut:

raw_data={'Alue': ['Helsinki', 'Turku', 'Tampere', 'Oulu'],
    'Myynti': [1321847, 852669, 1032199, 568230]}
 myynnit = pd.DataFrame(raw_data)
 myynnit.index = myynnit['Alue']
 myynnit

Yksinkertaisimmillaan luon vaakapylväskuvion myynnit-dataframen sisällöstä komennolla:

myynnit.plot.barh()

kuviot1

Kuvion muotoilu

Kuvion ulkoasuun vaikutan nopeimmin vaihtamalla tyyliä. Tarjolla olevien tyylien nimet näen komennolla plt.style.available. Tyylin vaihto esimerkiksi tyyliin ’ggplot’ onnistuu komennolla plt.style.use(’ggplot’).

Jos haluan räätälöidä kuviota yksilöidymmin, niin tallennan kuvion luonnin palauttaman Axes-objektin muuttujan arvoksi. Axes-objektin kautta voin muotoilla kuviota monin tavoin. Seuraavassa

  • otan käyttöön ’ggplot’-tyylin
  • järjestän myynnit suuruusjärjestykseen ennen plot.barh-komentoa
  • lisään kuvioon pääotsikon (title) ja poistan selitteen (legend) näkyviltä
  • tallennan kuvion luonnin palauttaman Axes-objektin ax-nimiseen muuttujaan
  • määritänn otsikon x-akselille (xlabel) ja tyhjennän y-akselin otsikon (ylabel)
  • muotoilen x-akselin luvut miljooniksi euroiksi (jos tähän liittyvä koodi vaikuttaa esoteerisiltä, niin siltä se minunkin mielestäni aluksi vaikutti).
plt.style.use('ggplot')
ax=myynnit.sort_values(by='Myynti').plot.barh
   (title='Myynti vuonna 2017', legend=False)
ax.set(xlabel='Miljoonaa euroa', ylabel='')
vals = ax.get_xticks()
ax.set_xticklabels(['{:.1f}'.format(x/1000000) for x in vals])

kuviot2

Seuraavaksi lisään dataframeen prosenttiosuudet kokonaismyynnistä ja esitän prosenttiosuudet pylväskuviona:

n=myynnit['Myynti'].sum()
myynnit['%'] = myynnit['Myynti'] / n * 100
ax=myynnit.sort_values(by='%')['%'].plot.barh
   (title='Osuus kokonaismyynnistä vuonna 2017', legend=False, 
   width=0.7, color='C0')
ax.set(xlabel='%', ylabel='')
vals = ax.get_xticks()
ax.set_xticklabels(['{:.0f} %'.format(x) for x in vals])

kuviot3

Lisämääritys width=0.7 vaikuttaa pylväiden leveyteen ja samalla pylväiden väliseen tyhjään tilaan. Oletus on width=0.5. Mitä suurempi arvo sitä lähempänä pylväät ovat toisiaan. Jos width=1, niin pylväät ovat kiinni toisissaan.

Outo värimääritys color=’C0′ vaatinee selitystä. Jos teen kuvion vain yhdestä dataframen sarakkeesta (esimerkiksi myynnit[’%’]) niin kyseessä ei ole enää dataframe, vaan series. Tällöin jokainen pylväs esitetään omalla värillään (vanhemmissa versioissa tällaista ominaisuutta ei ollut). Verkkokeskusteluiden perusteella on epäselvää, onko kyseessä virhe vai tarkoituksellinen ominaisuus. Esimerkin tapauksessa värit eivät mielestäni tuo lisäarvoa. Värimääritys color=’C0′ pakottaa pylväisiin käytössä olevan väripaletin ensimmäisen värin. Vastaavasti ’C1’ pakottaisi käytössä olevan väripaletin toisen värin jne.

Vierekkäset pylväät

Seuraavaksi luon dataframen, jossa on eri alueiden myynnit vuoden 2017 lisäksi myös vuodelta 2016 ja esitän molempien vuosien myynnit vierekkäisinä pylväinä:

raw_data={'Alue': ['Helsinki', 'Turku', 'Tampere', 'Oulu'], 
   'Myynti 2017': [1321847, 852669, 1032199, 568230],
   'Myynti 2016': [1203434, 923450, 1023563, 542399]}
myynnit2 = pd.DataFrame(raw_data)
myynnit2.index = myynnit2['Alue']
ax=myynnit2.plot.bar(figsize=(10,6), title='Myynti vuosina 
   2016 ja 2017', rot=0, width=0.7)
ax.set(xlabel='Alue', ylabel='Euroa')

kuviot4

Lisämääreenä määritin kuvion koon (figsize=(10,6)) ja käänsin luokka-akselin tekstit vaakaasentoon (rot=0).

Pinotut pylväät

Lisämääreellä stacked=True voin laatia pinotun pylväskaavion

ax=myynnit2.plot.bar(figsize=(10,6), title='Myynti 
   vuosina 2016 ja 2017', rot=0, width=0.7, legend=
   'reverse', stacked=True)
ax.set(xlabel='Alue', ylabel='Euroa')

kuviot5

Lisämääreellä legend=’reverse’ vaihdoin selitteen järjestyksen, jotta se vastaisi pinottujen pylväiden järjestystä.

Pinotut 100 % pylväät

Katsotaan vielä, miten laadin ristiintaulukoinnista 100 % pylväskuvion. Avaan datan, jossa on tietoja erään yrityksen työntekijöistä (muiden muassa sukupuoli ja koulutus) ja ristiintaulukoin koulutuksen ja sukupuolen. Kuvioksi sopii 100 % pinottu pylväskuvio:

df1 = pd.crosstab(df['koulutus'], df['sukupuoli'], margins = 
   True, normalize = 'columns')
df1.index = ['peruskoulu','2.aste','korkeakoulu',
   'ylempi korkeakoulu']
df1.columns = ['mies','nainen','yhteensä']

ax=df1[['mies','nainen']].transpose().plot.
   barh(stacked=True)
ax.set_xlabel('Prosenttia sukupuolesta')
ax.legend(loc='upper center', bbox_to_anchor=
   (0.5, 1.2), ncol=4)
vals = ax.get_xticks()
ax.set_xticklabels(['{:.0f} %'.format(x*100) 
   for x in vals])

kuviot6

Laadin ristiintaulukoinnin käyttämällä sukupuolta sarakemuuttujana (monilla on tapana ottaa selittävä muuttuja sarakemuuttujaksi). Tästä en kuitenkaan saa laadittua halutunlaista 100 % pylväskuviota ilman transpose()-toimintoa. Transpose() vaihtaa rivi ja sarakemuuttujat keskenään kuvion laadintaa varten.

Selitteen (legend) sijoittelukomento on taas melko esoteerisen näköinen. Pääset perille määreen bbox_to_anchor=(0.5, 1.2) tarkoituksesta kokeilemalla sulkeiden sisällä erilaisia lukuarvoja. Määre ncol=4 määrittää selitteen leveydeksi 4 saraketta (jokainen sarake sisältää yhden värin selitteen), joten värien selitykset sijoittuvat vierekkäin.

Figure-objekti

Plot-toiminnon myötä taustalla luodaan figure-objekti ja tämän sisälle axes-objekti. Voin tallentaa kuvion talletamalla figure-objektin:

plt.gcf().savefig('kuva.png')
plt.tight_layout()

gcf()-funktio palauttaa viimeksi luodun figure-objektin (get current figure). Tallennusformaatti määräytyy tiedostonimen tarkentimen (esim. png) perusteella. Matplotlib tukee yleisimpiä kuvaformaatteja.

plt.tight_layout() on tarpeen, jos tallennettu kuva ilman sitä näyttää reunoista liian ahtaalta.

Yhden figure-objektin sisälle on mahdollista luoda useita vierekkäin ja/tai allekkain sijoitettuja axes-objekteja. Yhden figure-objektin sisällä voi siis olla esimerkiksi useita pylväskuvioita.

Dataframe.plot

Tarkkaan ottaen hyödynsin tämän artikkelin kuvioiden luonnissa pandas-kirjaston dataframe.plot-toimintoa. Tämä toiminto kutsuu matplotlib.pyplot-kirjaston plot toimintoa. Lue lisää https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.plot.html

Lisää esimerkkejä pylväskuvioista löydät artikkelistani Pylväskaavioita Pythonilla.

Merkitsevyyden testaus Pythonilla

Tämän artikkelin ohjelmakoodin ja tulosteet löydät GitHubista:

https://github.com/taanila/tilastoapu/blob/master/p.ipynb

Jos kopioit koodia itsellesi, niin kannattaa käyttää GitHubia. Tästä artikkelista kopioidut koodit eivät välttämättä toimi oikein.

Oletan, että lukijalla on asennettuna Anaconda ja sen mukana tuleva Jupyter notebook.

Otoksessa havaitsemieni erojen ja riippuvuuksien tilastollista merkitsevyyttä voin arvioida laskemalla p-arvon. Pythonin scipy.stats-ohjelmakirjastosta löydän funktiot p-arvojen laskentaan.

Otan ensiksi käyttöön pandas ja scipy.stats -ohjelmakirjastot ja avaan esimerkkinä käyttämäni aineiston:

import pandas as pd
import scipy.stats as stats

df = pd.read_excel('http://taanila.fi/data1.xlsx', 
   sheet_name = 'Data')
df.head()

Korrelaatiokertoimen testaus

Iän ja palkan välisen pearsonin korrelaatiokertoimen ja siihen liittyvän 2-suuntaisen p-arvon saan funktiolla

stats.pearsonr(df['ikä'], df['palkka'])

Jos haluankin käyttää spearmanin järjestyskorrelaatiota, niin saan korrelaatiokertoimen ja 2-suuntaisen p-arvon funktiolla

stats.spearmanr(df['ikä'], df['palkka'])

Korrelaatiokertoimen testaamiseen liittyvistä funktioista löydät lisätietoa scipy.org -sivustolta:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html#scipy.stats.pearsonr

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.spearmanr.html

Ristiintaulukointi ja khiin neliö -testi

Esimerkiksi sukupuolen ja perhesuhteen väliseen ristiintaulukointiin liittyvän khiin neliö -testin testimuuttujan, p-arvon, vapausasteiden määrän ja odotetut frekvenssit saan funktiolla:

stats.chi2_contingency(pd.crosstab(df['sukup'], 
   df['perhe']))

Lisätietoa khiin neliö -testistä ja sen edeltävyysehdoista löydät scipy.org-sivustolta:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2_contingency.html#scipy.stats.chi2_contingency

Kahden riippumattoman otoksen t-testi

Jos haluan selvittää, onko miesten ja naisten palkkakeskiarvoissa eroa, niin erotan ensin miesten ja naisten palkat toisistaan (aineistossa 1=mies, 2=nainen)

a=df['palkka'][df['sukup']==1] #Mies
 b=df['palkka'][df['sukup']==2] #Nainen

Tämän jälkeen lasken t-testimuuttujan ja 2-suuntaisen p-arvon funktiolla

stats.ttest_ind(a, b, equal_var=False)

Yllä käytin erisuurten varianssien testiä (equal_var=False).

Lisätietoa riippumattomien otosten t-testistä ja sen edeltävyysehdoista löydät scipy.org-sivustolta:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html

Mann Whitney U -testi

Jos epäilen t-testin edeltävyysehtojen toteutumista, niin voin testata edellisen esimerkin Mann Whitney U-testillä:

stats.mannwhitneyu(a,b)

Tuloksena saan U-testimuuttujan ja p-arvon. Oletuksena saan 2-suuntaisen p-arvon puolikkaan. Lisätietoa löydät scipy.org-sivustolta:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mannwhitneyu.html

Yksisuuntainen varianssianalyysi

Jos haluan selvittää onko eri koulutuksen omaavien keskipalkoissa eroja, niin voin käyttää yksisuuntaista varianssianalyysiä (anova).  Ensiksi erotan eri koulutuksen omaavien palkat toisistaan:

k1=df['palkka'][df['koulutus']==1] #peruskoulu
k2=df['palkka'][df['koulutus']==2] #2. aste
k3=df['palkka'][df['koulutus']==3] #korkeakoulu
k4=df['palkka'][df['koulutus']==4] #ylempi korkeakoulu

Tämän jälkeen lasken anovan F-testimuuttujan ja p-arvon funktiolla:

stats.f_oneway(k1,k2,k3,k4)

Lisätietoa anovasta ja sen edeltävyysehdoista löydät scipy.org-sivustolta:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.f_oneway.html#scipy.stats.f_oneway

Kruskal-Wallis -testi

Jos epäilet varianssianalyysin edeltävyysehtojen  täyttymistä edellisessä esimerkissä, niin voit käyttää varianssianalyysin sijasta Kruskal-Wallis -testiä:

stats.kruskal(k1, k2, k3, k4)

Tuloksena saat H-testimuuttujan ja p-arvon.

Lisätietoa Kruskal-Wallis-testistä löydät scipy.org-sivustolta:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kruskal.html#scipy.stats.kruskal

Muita testejä

Lisää merkitsevyystestejä ja muita tilastollisia funktioita löydät scipy.org-sivustolta:

https://docs.scipy.org/doc/scipy/reference/stats.html