RACUNALNIŠ TVO O histogramih •is ■i' ■i' Primož Peterlin Histogram je grafični pripomoček za prikaz množice meritev zvezne številske spremenljivke. Prispevek predstavi konstrukcijo histograma, kumulativni histogram, razliko med histogramom in stolpičnim diagramom, težave s histogrami, za konec pa se pomudi še ob ocenjevanju porazdelitve z jedri. Od podatkov do histograma Pri fizikalnih meritvah pogosto merimo vrednost ene količine v odvisnosti od druge, npr. napetost na kondenzatorju v odvisnosti od časa. Še preprostejše pa so meritve, pri kateri merimo vrednosti ene same količine, npr. telesno višino učenčev v razredu ali število jedrskih razpadov, zaznanih z Geiger-Mullerjevo cevjo. V statistiki takšnim primerom, pri katerih imamo opravka z eno samo spremenljivko (telesna višina, število razpadov v časovni enoti), pravimo univariatni, za razliko od bivariatnih, kadar sta spremenljivki dve, ali v splošnem multivariatnih. Vsi učenči niso enako visoki in v izbranem časovnem intervalu ne razpade vedno enako število jeder, zato tako zbrani podatki navadno niso vsi enaki. Za zgled smo zbrali rezultate 49-ih učencev dveh para-lelk 6. razreda pri skoku v daljavo z mesta; rezultati so podani v čentimetrih. 135 168 160 166 130 171 148 152 156 120 176 139 189 115 130 135 140 134 180 180 125 106 141 169 193 129 139 130 165 149 120 148 140 95 150 176 184 159 152 169 185 147 150 190 175 120 149 155 174 Množica podatkov je nepregledna, zato želimo informacijo strniti in na preprost nacin predstaviti nje- ne glavne značilnosti. Dva uporabna parametra za opis množice podatkov sta povprečje in standardni odklon. Prvi pove »težišče« podatkov, drugi pa, koliko se podatki med seboj razlikujejo. Standardni odklon je uporabno merilo za opis raznolikosti podatkov, včasih pa želimo vedeti še kaj več o tem, kako so podatki porazdeljeni. Ali obstaja ena takšna vrednost, okrog katere so izmerjene vrednosti posejane pogosteje kot sičer, ali pa morda dve ali čelo več takih vrednosti? So vrednosti okoli povprečja posejane simetrično ali niso? Pripomoček, s katerim lahko dobimo približno grafično sliko o porazdelitvi meritev univariatne številske spremenljivke, je histogram. Kako se lotimo priprave histograma? Podatke razvrstimo v razrede ali predalčke. V zgornjem zgledu, denimo, zberemo skupaj rezultate med 91 in 100 čm, potem med 101 in 110 čm itd. Zgornjo in spodnjo mejo postavimo tako, da zajamemo vse podatke, širino predalčka pa izberemo tako, da se kar najbolj pokaže oblika porazdelitve. Prevelika širina predalčka bo morda zgladila in skrila kakšno sičer morda zanimivo lastnost porazdelitve, ob premajhni pa bo v posameznem predalčku premalo podatkov in prevladale bodo naključne razlike. Ko bomo končali, bo število v posameznem predalčku govorilo o tem, kako pogosto se določene vrednosti pojavljajo v vzorču; pogostosti se s tujko pravi tudi frekvenca (ki pa s frekvenčo pri nihanju in valovanju nima neposredne zveze). Razvrščanje podatkov v predalčke je na srečo posel, ki ga dobro obvladajo programi za obdelavo podatkov. Za zglede v prispevku bomo uporabili prosti programski paket GNU R (http://www. r-project. org/), o katerem smo v Preseku že pisali [2]. # datoteka z dolžinami skokov x <- scan("skoki.txt") hist(x, main=NULL, xlab = "dol\u017Eina (cm)", ylab = "frekvenca") PRESEK 42 (2014/2015) 6 23 RACUNALNIŠ TVO —^ Programček je tako kratek, da o njem skoraj ni kaj povedati. Z ukazom scan() v vektor x preberemo podatke iz datoteke skoki .txt, kjer je v zapisana po ena vrednost v vrstici. Klic funkcije hist(), ki ji vektor x podamo kot argument, opravi vse ostalo: vrednosti razvrsti v razrede in izriše histogram. Ker imamo podpis pod sliko, se odrečemo naslovu histograma (main=NULL), z izbirama xlab in ylab pa nastavimo oznaki na abscisi in ordinati. V oznakah lahko uporabimo tudi znake izven nabora ASCII; pri-klicemo jih s kodami ISO10646/Unicode. Rezultat vidimo na sliki 1a. Vidimo, daje hist() samodejno izbral povsem razumne meje intervalov. Če z njegovo izbiro ne bi bili zadovoljni, bi lahko z izbiro breaks= sami dolocili meje razredov. Prikaz frekvence (števila v posameznem podatkovnem razredu) ni edina mogoČca predstavitev histograma. Če število v vsakem razredu delimo s številom vseh meritev v vzorcu, dobimo relativne frekvence, npr. 1/49 = 0,02 ipd. Če pa relativno frekvenco delimo še s širino razreda, tak histogram prikazuje gostoto verjetnosti. Programcek v tem primeru popravimo tako, da funkciji hist() dodamo izbiro probabi l i ty=TRUE: 1 2 3 4 5 6 # datoteka z dolžinami skokov x <- scan("skoki.txt") hi st(x, probabili ty=TRUE, main=NULL, xlab = "dol\u017Eina (cm)", ylab = "gostota verjetnosti") Starost Frekvenca 0-4 28 5-9 46 10-15 58 16 20 17 31 18-19 64 20-24 149 25-59 316 60+ 103 Rezultat je prikazan na sliki 1b. Gostota verjetnosti je definirana tako, da je skupna plošcina vseh stolpcev natanko 1. TABELA 1. Udeleženci prometnih nesrec v londonskem okrožju Harrow v letu 1985. Razlika med obema je morda najbolj očitna v primeru, ko razredi niso enako široki (zgled je izposojen iz učbenika statistike, [2, str. 25]). V londonskem okrožju Harrow so za leto 1985 zbrali statistiko udeležencev prometnih nesreč po starosti (tabela 1). Podatke najprej prikažemo s stolpičnim diagramom. 1 vrednosti <- c(28, 46, 58, 20, 31, 2 64, 149, 316, 103) 3 barplot(vrednosti, 4 names.arg = c("0-4", "5-9", 5 "10-15", "16", "17", 6 "18-19", "20-24", "25-59", 7 "60+"), 8 xlab = "Starost (leta)", 9 ylab = "\u0160tevilo") Histogram in stolpični diagram Histogram pogosto zamenjujejo s stolpicnim diagramom. Kljub temu, da so pri obeh podatki predstavljeni s stolpci, pa je razlika med njima precejšnja: StoplpiČcni diagram prikazuje frekvence na diskretni osi kategorialne spremenljivke, bodisi nominalne (npr. moški/ženske) bodisi ordinalne (npr. stopnja izobrazbe). ■ Histogram je približek porazdelitve po zvezni spremenljivki. Programcek je spet zelo preprost. Vektor vrednosti vsebuje frekvence (desni stolpec v tabeli 2), stolpicni diagram pa izrišemo s funkcijo barplot(), ki ji podamo ta vektor. Argument names.arg je vektor z oznakami za posamezne stolpce, xl ab in yl ab pa oznaki osi. Rezultat programa je stolpicni diagram, prikazan na sliki 2. Diagram pravzaprav ne pove vec od tabele 1 in ne odraža porazdelitve udeležencev prometnih nesrec po starosti. Vidimo lahko, denimo, da je število odraslih udeležencev prometnih nesrec (torej tistih v starostni skupini 25-59 let) desetkrat 24 PRESEK 42 (2014/2015) 6 24 RACUNALNIŠ TVO □ I I I I I l 00 140 180 dolžina (cm) (a) o s o C o CU o > rt Ln 0 <3 o" 0 0 <3 <3 □ I I I I I I 00 140 180 dolžina (cm) (b) SLIKA 1. Histogram porazdelitve skokov v dolžino. tolikšno kot število udeležencev, starih 17 let. Vendar prva skupina zajema dosti vecji delež populacije kot druga, zato nam ta podatek sam po sebi ne pove veliko. Predstavimo zdaj podatke iz tabele 2 kot histogram. Za razliko od primera na sliki 1 so predalčki tu različno široki. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 meje <- c(0, 5, 10, 16, 17, 18, 20, 25, 60, 80) vrednosti <- c(28, 46, 58, 20, 31, 64, 149, 316, 103) mojhist <- list(breaks=meje, counts=vrednosti, densi ty=vrednosti/di ff(meje), xname=NULL) class(mojhist) <- "histogram" plot(mojhist, xlab="Starost (leta)", ylab="\u0160tevilo na leto starosti", mai n=NULL) V tem primeru ne potrebujemo klica funkcije hist(), ki razvrsti podatke po predalčkih, saj je nekdo to že opravil namesto nas. V vektorju mej e so shranjene meje razredov (najvišji interval smo omejili na 80 let), v vektorju vrednosti pa frekvence posameznih razredov. Zatem sestavimo seznam mojhist z elementi breaks, kateremu podamo meje; counts, kateremu podamo frekvence; density, kateremu podamo vrednosti, deljene s širino razreda; (diff(meje) vrne vektor razlik med zaporednimi elementi vektorja meje, kar so ravno širine razredov); naslova (xname) pa ne nastavimo. Potem uporabimo predmetno naravo jezika R in seznamu mojhi st, ki smo ga ravnokar ustvarili takega, da že ima pravilno strukturo histograma, priredimo razred histogram. Končno histogram mojhist izrišemo s funkcijo plot(). Izbire main, xlab in ylab imajo enak pomen kot prej. Rezultat je na sliki 3. Diagram se precej razlikuje od tistega na sliki 2. Vidimo lahko, da so mladostniki približno trikrat pogosteje udeleženci prometnih ne-srec kot odrasli (ali tudi kot otroci). Slika 3 tudi nazorno pokaže, da je na histogramu frekvenci sorazmerna ploščina posameznega stolpca, ne pa njegova višina. Pri razredih enake širine sta plošcina in višina stolpca resda premo sorazmerni, kar lahko zavede. Kumulativni histogram Iz histograma na sliki 1 lahko neposredno preberemo podatek o frekvenci posameznega razreda, denimo, koliko ucencev je skocilo med 140 in 149 cm. Vcasih pa nas zanima drugacno vprašanje: denimo, koliko ucencev je skocilo manj kot 150 cm. Odgovor lahko seveda izracunamo, tako da seštejemo frekvence v razredih 90-99 cm, 100-109 cm in tako dalje do 140-149 cm. Kumulativni histogram (slika 4b) pa omogoca, da tak podatek preberemo neposredno iz diagrama. Poseben primer obrnjenega kumulativnega histograma so tudi krivulje preživetja, ki se uporabljajo v biomedicinskih vedah [4]. 1 x <- scan("skoki.txt") 2 h <- hi st(x, plot=FALSE) 3 h$counts <- cumsum(h$counts) 4 plot(h, main=NULL, 5 xlab = "dol\u017Eina (cm)", 6 ylab = "frekvenca") PRESEK 42 (2014/2015) 6 21 RACUNALNIŠ TVO o o o o o o SLIKA 2. Stolpicni diagram števila udeležencev prometnih nesrec po starosti. 0-4 5-9 10-15 16 1 7 1 8-19 20-24 25-59 60+ Starost (leta) SLIKA 3. Histogram števila udeležencev prometnih nesreč po starosti. 20 -1- 40 Starost (leta) 60 80 Oglejmo si še, kako v R pripravimo kumulativni histogram. Enako kot prej v vektor x preberemo podatke iz datoteke. Dve razliki pa sta pri funkciji hist(). Prva je ta, da smo ji podali izbiro plot=FALSE, s katero smo zahtevali, da histogram sicer izračuna (določi podatkovne razrede in vanje razvrsti podatke), a ga ne izriše. Druga pa je ta, da smo rezultat izračuna histograma shranili v spremenljivko h. Iz zgleda z razredi neenake širine že vemo, da je ta spremenljivka seznam; skladnja h$counts nam vrne element seznama counts, ki je vektor s frekvencami. V tretji vrstici histogram pre- tvorimo v kumulativni histogram s klicem funkcije cumsum(). Ce ji podamo vektor dolžine n, bo vrnila vektor iste dolžine s kumulativnimi vsotami: na prvem mestu bo kar prvi element podanega vektorja na drugem mestu vsota prvih dveh elementov podanega vektorja, in tako dalje do zadnjega elementa, kjer bo vsota vseh elementov podanega vektorja. Vektor s kumulativnimi vsotami shranimo kot element h$counts histograma. Tako spremenjen histogram zdaj narišemo z ukazom plot(), ki mu podamo iz-racunani histogram h, izbire main, xlab in ylab pa imajo že znani pomen. 26 PRESEK 42 (2014/2015) 6 26 RACUNALNIŠ TVO I I I I I I 100 140 180 dolžina (cm) (a) 0 4 i i i i i i 100 140 180 dolžina (cm) (b) SLIKA 4. Histogram (a) in kumulativni histogram (b) porazdelitve skokov v dolžino. Težave s histogrami V prvem zgledu smo nekoliko zlahka odpravili izbiro intervala in določitev števila razredov. Cas je, da priznamo, da sta prav ti dve izbiri srž težav s histogrami. Posebej v primeru, če podatkov ni veliko, je histogram odvisen od izbire teh dveh parametrov. Oglejmo si najprej zgled, kako na histogram vpliva izhodišče razredov. Levi diagram na sliki 5 že poznamo, pri desnem (slika 5b) pa smo izbrali razrede tako, da zaobjamejo vrednosti 95-104, 105-114 itd. V programčku smo to izvedli z izbiro breaks=, pri kateri smo uporabili funkcijo seq(). Tej smo podali tri parametre: prvi element, zadnji element in korak; funkcija vrne vektor z zaporedjem, ki sledi podanemu pravilu. 1 x <- scan("skoki.txt") 2 hist(x, main=NULL, 3 breaks = seq(95, 205, 10), 4 xlab = "dol\u017Eina (cm)", 5 ylab = "frekvenca") Vnaprej lahko uganemo, da večja širina razreda zgladi histogram, kar lahko vidimo tudi na sliki 6. Levi histogram uporablja privzete meje razredov (breaks = seq(90, 200, 10)), desni pa dvakrat tolikšno širino razredov (breaks = seq(90, 210, 20)). Problem izbire optimalnega števila razredov so precej preučevali in različni raziskovalci so prišli do različnih formul za optimalno število razredov. Altman kot praktični nasvet navaja [1], da je navadno dovolj 8-15 razredov, razen če je podatkov zelo veliko. Med bolj znanimi so še Sturgesova formula, k = riog2 n + 11, kjer je n število podatkov, k število razredov, \x1 pa označuje zaokroževanje navzgor, in Scottova formula h = 3,5a/n1/3, kjer je a standardni odklon vzorca, h pa širina razreda. Kateri od histogramov je pravi? Pravega ali najboljšega histograma ni. Ne poznamo postopka, s katerim bi za poljubno porazdelitev vhodnih podatkov izračunali najboljši histogram. Na nas je, da s po-skušanjem in spreminjanjem izhodišča in širine razredov izračunamo histogram, ki je sprejemljiv. Na srečo si lahko pomagamo z računalnikom, kar vsakokratno razvrščanje podatkov v razrede napravi skoraj hipno. Zato je histogram kljub naštetim pomanjkljivostim še vedno uporabno orodje za kvalitativno očeno eksperimentalnih porazdelitev. Ocenjevanje gostote verjetnosti z jedri Histogram ima zaradi svoje enostavne konstrukčije in interpretačije zagotovo svoje prednosti, vseeno pa se moramo glede na vse prej omenjene težave s histogrami vprašati, ali ne obstaja kakšen boljši način za očeno porazdelitve v vzorču, pridobljenem s poskusom. Obstaja. Metoda je poznana kot očenjeva-nje gostote z jedri (angl. kernel density estimation). Matematično je znatno bolj zapletena in preobsežna za ta članek. Osnovna zamisel pa je preprosta in jo bomo nakazali. Za zgled si oglejmo vzoreč 12-ih meritev (vrednosti nalašč ni preveč, da je primer preglednejši): 2,064 2,212 2,351 2,409 2,459 2,639 2,656 2,673 3,350 3,373 3,599 3,861 PRESEK 42 (2014/2015) 6 27 RACUNALNIŠ TVO ra o O (a) (b) (c) SLIKA 7. Ocenjevanje gostote verjetnosti z Gaussovimi jedri; (a) premajhno glajenje, (b) preveliko glajenje, (c) optimalno glajenje. Rdece črtkane krivulje podajajo prispevke posameznih jeder, crna neprekinjena crta pa na osnovi jeder dobljeno gostoto verjetnosti. Črtice na notranji strani abscisne osi označujejo vrednosti meritev. V zgledu prispevek vsake posamične meritve opi-šimo z Gaussovo porazdelitvijo, ki ima vrh pri vrednosti te meritve. Matematični funkciji, s katero opišemo prispevek posamične vrednosti, pravimo jedro. Skupno porazdelitev potem dobimo kot seštevek posamičnih prispevkov. Porazdelitev je odvisna tudi od širine Gaussovega jedra, ki jo podaja standardni odklon c. Na sliki 7 so predstavljene porazdelitve gostote verjetnosti, dobljene z Gaussovimi jedri z različnimi vrednostmi c: c = 0,1 (a), c = 0,5 (b) in c = 0,25 (č). V R očeno gostote z jedri izračuna fukčija densi ty(). Na sliki 7 vidimo, da širina jedra močno vpliva na očeno gostote verjetnosti. To še ni vse: brez obrazložitve smo za jedro vzeli Gaussovo funkčijo, lahko pa bi tudi kakšno drugo. Smo torej sploh kaj na boljšem kot pri histogramih, ki smo jim očitali preveliko subjektivnost? Malo bolje je vseeno. Načeloma je naloga preprosta: optimalna, z jedri očenjena porazdelitev je tista, ki se čimbolj ujema s pravo; težava pa je v tem, da prave porazdelitve ne poznamo, ampak bi jo radi šele ugotovili. Kljub vsemu obstajajo postopki, ki v asimptotičnem približku vodijo k optimalni jedrni funkčiji in optimalni širini. Ne najnovejši, pa še vedno prečej bran učbenik s tega področja je [3], področje pa se še vedno razvija. Za koneč povzamimo dobre in slabe lastnosti obeh pristopov. Histogram je preprosto konstruirati s svinčnikom in milimetrskim papirjem, preprosto ga je interpretirati in kljub težavam z arbitrarnostjo izbire izhodišča in širine razredov večinoma nudi grobo očeno za porazdelitev izmerjenih vrednosti. Oče-na gostote z jedri je matematično bolj kompleksna, računsko zahtevnejša, nudi pa nekoliko boljšo očeno porazdelitve. Dostopnost zmogljivih računalnikov ter izvedbe v večini programskih jezikov in paketov pomenita, da je ta metoda, nekoč omejena na raziskovalne laboratorije, dostopna vsakomur. Zato je dobro, da tudi razumemo, kako deluje. Literatura [1] D. G. Altman, Practical statistics for medical research, London: Chapman & Hall, 1991. [2] P. Peterlin, Obdelava meritev in risanje grafov z R, Presek 37, 24-30, 2010. [3] B. W. Silverman, Density estimation for statistics and data analysis, London: Chapman & Hall, 1986. [4] J. Stare, Krivulje preživetja, Medicinski razgledi 40,173-181, 2001. _XXX PRESEK 42 (2014/2015) 6 29