Verohallinnon avoimen datan tarkastelua

Tässä tekstissä tutkin Verohallinnon avointa dataa, josta löytyvät tiedot julkisten yhteisöjen tuloverotuksesta vuodelta 2014. Datasta löytyvät mm. y-tunnus, verotuskunta, verotettavat tulot, maksuunpannut verot, ennakkoverot, veronpalautukset sekä jäännösverot. Hyödynnän lisäksi verovelkarekisteriä, josta löytyy tieto, onko yrityksellä verovelkaa yli 10000 euroa. Haluan erityisesti selvittää, onko jäännösverolla ja verovelalla yhteyttä, eli toimiiko jäännösvero indikaattorina tulevaisuuden verovelalle. Suuri jäännösvero voi kertoa yrityksen maksuvaikeuksista tai siitä, että yrityksen tulos on ylittänyt odotukset, jolloin ennakkoveron olisi pitänyt olla suurempi. Toisaalta tällöin yritys olisi voinut maksaa puuttuvat ennakkoveronsa ennakontäydennysmaksuna, jolloin korollista jäännösveroa ei olisi syntynyt. Lisätietoa ennakkoverotuksesta löytyy mm. Verohallinnon sivuilta täältä.

Rajataan tarkastelu aluksi yrityksiin, joilla on verotettavaa tuloa yli 10000 euroa. Näin saadaan aineistoa rajattua suurempiin yrityksiin ja niihin, joilla on jäännösveroa. Lasketaan sitten yhteisö- ja jäännösveroprosentti seuraavasti:

R-koodi:

data <- read.csv(file.choose(), header=F,sep=";",dec=",",skip=1)
colnames(data) <- c("Verovuosi","Ytunnus","Verovelvollisennimi","Verotuskunta","Verotettavatulo","Maksuunpannutverotyhteensa","Ennakotyhteensa", "Veronpalautus","Jaannosvero")

data1 <- data[data$Verotettavatulo>10000&data$Maksuunpannutverotyhteensä>0,]

data1$verotusuhdekaikki <- data1$Maksuunpannutverotyhteensa/data1$Verotettavatulo*100

data1$jaannossuhdekaikki <- data1$Jaannosvero/data1$Maksuunpannutverotyhteensä*100

Katsotaan yhteisöjen tuloveroprosenttien persentiilejä ja nähdään, että mukana on myös yrityksiä, joilla yhteisöveroprosentti on reippaasti yli 20.

Katsotaan vielä histogrammia. Rajataan tarkastelu yrityksiin, joilla tuloveroprosentit ovat välillä 15 ja 30. Yritysten määrä tippuu 57816:sta 57113:een.

Lasketaan seuraavaksi kunnille prosenttiosuudet, joista selviävät, kuinka monta prosenttia kunnan yrityksistä on saanut yli 50%:n verojäämän. Tätä varten tehtään funktio, jota käytetään yhdessä lapply-komennon kanssa. Ennen funktiota muodostetaan uusi dataframe, jossa ovat vain kuntien nimet. Tätä varten pitää Verotuskunta-muuttuja ensin jakaa kahteen osaan, koska se sisältää sekä kuntakoodin että nimen. Tällainen kahtiajako on kätevä tehdä reshape-pakkauksen transform-komennolla.

install.packages("reshape")  # Tehdään muuttujan kahtiajako

library(reshape)

data3 = transform(data3, verotuskunta = colsplit(data3$Verotuskunta, split = " ", names = c('numeric', 'nimi','koodi')))

data4 <- data.frame(unique(data3$verotuskunta.nimi)) # Muodostetaan uusi dataframe, jossa kuntien nimet vain kerran

colnames(data4) <- "Kunta"

data4$Kunta <- as.character(data4$Kunta)

osuus <- function(i) {                                   # Tehdään funktio

a <- length(data3$verotuskunta.nimi[data3$verotuskunta.nimi==i])

b <- length(data3$jaannossuhdekaikki[data3$verotuskunta.nimi==i&data3$jaannossuhdekaikki>=50])

c=b/a*100

return(c)

}

data4$Osuus <- lapply(data4$Kunta,osuus) # Lasketaan prosenttiosuudet joka kunnalle

Tehdään karttakuva, josta lasketut prosenttiosuudet selviävät. Kuvaa varten täytyy ensin hakea koordinaatit, jotka löytyvät esim. täältä. Sitten datat yhdistetään kunnan perusteella. 

R-koodi:

koordinaatit<- read.csv(file.choose(),header=TRUE,dec=".",sep=";")

colnames(koordinaatit)[2] <- "Kunta"

datayhd <- merge(koordinaatit,data4, by="Kunta")

datay <- datayhd[!duplicated(datayhd$Kunta),]

datay$Osuus <- as.numeric(datay$Osuus)

g <- sample(nrow(datay),0.05*nrow(datay))    #Otetaan kuvaa varten 5% havainnoista

datas <- datay[g,]

 map <- get_map(location = 'Finland', zoom = 6)  #Ladataan Suomen kartta 

  mapPoints <- ggmap(map)+

   geom_point(aes(x=Longitude,y=Latitude,size=Osuus),data=datas,alpha=0.5,colour="blue")

  mapPoints

mapPointsLegend <-mapPoints +scale_size(range=c(7,25)) +ggtitle("Jäännösveroprosentti >50: Osuudet Suomen kunnissa")

mapPointsLegend

Katsotaan vielä muutamia suurimpien ja pienimpien osuuksien kuntia. Helsingissä Osuus on ~25%. 

head(datay[order(datay$Osuus),])

          Kunta Postal.Code Latitude Longitude Osuus

123    Hailuoto       90480  65.0000   24.7167     0

624   Kivijärvi       43801  63.0667   25.0500     0

705    Kumlinge       22820  60.2667   20.7833     0

896     Luhanka       19950  61.7833   25.7000     0

899  Lumparland       22630  60.1167   20.2500     0

1298  Savukoski       98800  67.2925   28.1581     0

head(datay[order(-datay$Osuus),])

             Kunta Postal.Code Latitude Longitude     Osuus

909      Merijärvi       86220  64.3000   24.4500 100.00000

118           Geta       22340  60.3833   19.8500  75.00000

1205    Rautavaara       73901  63.4833   28.3000  50.00000

1047 Pelkosenniemi       98500  67.1108   27.5106  41.66667

117          Föglö       22710  60.0167   20.4167  37.50000

1536         Vaala       91700  64.4333   26.8000  36.84211

Seuraavaksi tutkitaan, onko suuren jäännösveroprosentin yrityksillä suurempi todennäköisyys saada verovelkaa. Rajaan tutkittavan joukon verotettavan tulon mukaan välille 50000-80000. Tästä joukosta valitsen ne, joilla on yli 99%:n jäännösvero ja ne, joilla on alle 1%:n jäännösvero. Näistä kummastakin joukosta otan sitten 30 kappaleen satunnaisotannan. 

data33 <- data3[data3$Verotettavatulo>=50000 & data3$Verotettavatulo<=80000,]

datajaannos99 <- data33[data33$jaannossuhdekaikki>99,]

datajaannos1 <- data33[data33$jaannossuhdekaikki <1,]

a <- sample(nrow(datajaannos99),30)

b <- sample(nrow(datajaannos1),30)

datajaannos99t <- datajaannos99[a,]

datajaannos1t <- datajaannos1[b,]

Yli 99%:n jäännösveron otoksessa verovelkaisia oli 7 kappaletta, eli 23.3%:a kaikista. 

Alle 1%:n jäännösveron otoksessa verovelkaisia ei ollut yhtään kappaletta. 

Tehdään vielä Fisherin exact testi. Testillä voidaan testata, onko populaatioiden dikotomisten muuttujien arvojen suhteellisissa osuuksissa eroja.

data = matrix(c(7,23,0,30), ncol=2)
fisher.test(data,alternative = "two.sided")

Fisher's Exact Test for Count Data
data:  data
p-value = 0.01054
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.649705      Inf
sample estimates:
odds ratio 
 Inf 

Ratkaistaan myös luottamusväli yli 99%:n ryhmän verovelallisten osuudelle.

prop.test(x=7,n=30,conf.level = 0.95)

1-sample proportions test with continuity correction

data:  7 out of 30, null probability 0.5

X-squared = 7.5, df = 1, p-value = 0.00617

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

 0.1063502 0.4270023

sample estimates:

        p 

0.2333333 

Laskin verovelallisten osuuden lopuksi myös 45-55 %:n jäännösveroryhmälle. 30:n yrityksen otoksessa oli verovelallisia 3 kappaletta, eli 10%:a kaikista.

Verovelkarekisterissä näkyvät tiedot vain viimeiseltä 6 kuukaudelta, joten käytetyn datan (2014) ja verovelkadatan välillä on yli vuoden ero. Tämä antaa viitteitä siitä, että suuri jäännösveroprosentti tänään voi merkitä suurempaa verovelkaisuuden todennäköisyyttä tulevaisuudessa. Ei ole tietenkään suuri yllätys, että alle 1%:n jäännösverolaisilla ei ollut verovelkaa, koska hehän ovat jo maksaneet veronsa ennakkoveroina, jolloin verovelkaa ei voi sillä hetkellä syntyä. Tämä tulos on kuitenkin hyödyllinen, koska se ennakoi verovelvollisen käytöstä tulevaisuudessa. 

R ja suuret aineistot

Esittelen nyt suurten aineistojen käsittelyä R:llä. Käytän harjoitusaineistona Englannin ja Walesin asuntojen hintadataa, joka löytyy täältä. Datassa on tietoja asuntojen hinnoista Englannissa ja Walesissa vuodesta 1995 lähtien. Koko data sisältää yli 22 miljoonaa riviä ja 16 muuttujaa. Muuttujina löytyvät mm. asunnon hinta, uusi/vanha asunto, tontti oma/vuokralla, myyntipäivämäärä ja asuntotyyppi. Asuntojen koosta ei ole kuitenkaan dataa. Olen kiinnostunut tutkimaan, onko kalleimpien asuinpaikkojen keskihinnoissa ollut eroja vuosina 2014 ja 2015. Rajaan tarkastelun omakotitaloihin, joilla on oma tontti.

Aloitetaan lataamalla data R:ään. Tähän kannattaa käyttää data.table-pakettia, josta löytyy fread-komento. Tämä on huomattavasti nopeampi tapa ladata suuri määrä dataa kuin esim. read.csv. Tarvitsen suuresta datamäärästä dataa vuosille 2014 ja 2015, joten jätän osan datasta lataamatta käyttämällä skip-komentoa. Dataa latautuu 2.5 miljoonaa riviä. Tällaisen määrän dataa R pyörittää vielä hyvin 6GB:n keskusmuistilla.

library(data.table)
data <- fread(file.choose(),header=F,skip=19000000)
colnames(data) <- c("ID","Price","Date","Postcode","Propertytype","Oldnew","Duration","PAON","SAON","Street","Locality","Towncity","District","County","PPD","RS")
data$Price <- as.numeric(data$Price)
data$Date <- as.Date(data$Date, format="%Y-%m-%d")
data$Year <- as.numeric(format(data$Date,"%Y")) # Teen vuosi-muuttujan, jotta vuodet on helppo valita jatkossa.

Nyt kun dataan on tehty sopivat muunnokset, lasketaan ensin uusien omakotitalojen keskihinnat vuosille 2014 ja 2015.

mean(data$Price[data$Year=="2014"  & data$Oldnew=="Y"&data$Duration=="F"&data$Propertytype=="D"]) #Keskihinnat 2014 ja 2015
mean(data$Price[data$Year=="2014"  & data$Oldnew=="Y"&data$Duration=="F"&data$Propertytype=="D"])

Keskihinnat ovat 2014: ~358000 ja 2015: ~373000

Tarkastellaan seuraavaksi keskihintoja maakuntatasolla.

data2015<- do.call(data.frame,aggregate(data$Price[data$Year==2015&data$Propertytype=="D"&data$Duration=="F"],list(data$County[data$Year==2015&data$Propertytype=="D"&data$Duration=="F"],data$Oldnew[data$Year==2015&data$Propertytype=="D"&data$Duration=="F"]),function(x) c(mean=mean(x),length=length(x))))

data2014<- do.call(data.frame,aggregate(data$Price[data$Year==2014&data$Propertytype=="D"&data$Duration=="F"],list(data$County[data$Year==2014&data$Propertytype=="D"&data$Duration=="F"],data$Oldnew[data$Year==2014&data$Propertytype=="D"&data$Duration=="F"]),function(x) c(mean=mean(x),length=length(x))))

Otetaan tarkastelusta pois alueet, joista on alle 30 havaintoa.

data2015l <- subset(data2015, data2015$Length>=30,)
data2014l <- subset(data2014, data2014$Length>=30,)

Järjestetään data keskiarvon mukaan suurimmasta pienimpään.

data2015l <- data2015l[order(-data2015l$Mean),]
data2014l <- data2014l[order(-data2014l$Mean),]

Katsotaan heatmappia, jossa on alhaisen ja korkean keskihinnan maakuntia vuodelta 2015.

datax <- data2015l[c(1:6),]
datay <- data2015l[c(184,189,185,200,191,195.192),]
data5 <- rbind(datax,datay)

library(ggplot2)
ggplot(data = data5, aes(x = Oldnew, y=County, fill=mean)) + 
geom_raster() + scale_fill_gradient2(low="red", mid="yellow", high="green",na.value="grey") +
labs(fill='Keskihinta',title="Asuntojen keskihinnat")

Tehdään seuraavaksi  kuva, jossa näkyvät neljä keskihinnaltaan kalleinta uusien omakotitalojen aluetta vuodelta 2015.

data2015 <- subset(data2015l,data2015l$Oldnew=="Y",)
data2015r <- data2015[,"Mean",drop=FALSE]  #drop=FALSE pitää data.frame -muodon,
row.names(data2015r) <- data2015$County  
data2015tm <- t(as.matrix(data2015r[1:4,,drop=FALSE]))


b <- barplot(data2015tm,,col=c("3"),las=1,yaxt="n",ylab="£",cex.names=0.62,,main="Englannin ja Walesin kalleimmat asuinalueet vuonna 2015 \n (Vertailu: Uusien omakotitalojen keskihinnat)")
axis(2,seq(0,1600000,by=200000))
legend("topright",legend=c("Keskihinta","Kaikkien uusien omakotitalojen keskihintataso"),col=c("3","1"),pch=c(19,NA),lty=c(NA,1))
text(b,60000,round(data2015tm,0))
abline(h=372630,lwd=2)

Katsotaan lopuksi, miten näiden alueiden keskihinnat ovat kehittyneet aiemmasta vuodesta.

data20157<- data2015[1:4,c("County","Mean")]
data2014 <- subset(data2014l,data2014l$Oldnew=="Y",)
yhd <- merge(data20157,data2014,by="County",all.y=F)
yhd$erotus20152014 <- yhd$Mean.x -yhd$Mean.y
yhd$muutosprosentti <- (yhd$erotus20152014/yhd$Mean.y)*100
colnames(yhd) <-c("County","mean2015","Oldnew","mean2014","length","erotus20142014","muutosprosentti")
yy <- yhd[,"muutosprosentti",drop=FALSE]
row.names(yy) <- yhd$County

muutosprosentti <- t(as.matrix(yy))
b <- barplot(muutosprosentti,col=c("5"),las=1,yaxt="n",ylab="%",cex.names=0.62,main="Keskihintojen muutos vuosina 2014-2015")
axis(2,seq(-20,25,by=1))
legend("topright",legend=c("Muutosprosentti"),col=5,pch=19)

text(b,1,round(muutosprosentti,2))

Vaihtelua vuosien välillä on siis melkoisen paljon. Ongelmana tässä datassa on, ettemme tiedä asuntojen kokoa. Asuntojen hintojen kehitystä olisi luontevinta tutkia vertailemalla neliöhintoja alueittain. Suuret hintavaihtelut selittyvät luultavasti osittain sillä, että tutkimme uusia omakotitaloja. Tällaisia taloja myytiin esim. Suur-Lontoossa 223 kappaletta vuonna 2015. Windsorissa vastaava luku oli vain 35. Tällaisilla kalliilla ja harvinaisilla taloilla hintavaihtelut ovat luultavasti suurempia kuin edullisemmilla taloilla, joita on markkinoilla tarjolla koko ajan. Alla on taulukko, jossa on muutosprosentteja vanhoille omakotitaloille kalleimmilla alueilla. Näemme, että luvut ovat Suur-Lontoolle, Poolelle, Surreylle ja Windsorille huomattavasti maltillisemmat. Lisäksi Count-sarakkeiden myydyistä kappalemääristä näemme, että Suur-Lontoossa ja Windsorissa myydään runsaasti enemmän vanhoja kuin uusia omakotitaloja, mikä ei tietenkään ole yllätys. Windsorista käytetyn omakotitalon sai vuonna 2015 halvimmillaan 94000 punnan hintaan, kun taas uudesta joutui maksamaan 325000 puntaa. Kallein myyty omakotitalo maksoi siellä 6.4 miljoonaa puntaa.

2015-2014 muutosprosentit vanhoille omakotitaloille

County	Mean2015	Count2015	Mean2014	Count2014	Muutosprosentti
GREATER LONDON	947019	5418	917501	5651	3.22
WINDSOR AND MAIDENHEAD	851027	607	812147	640	4.79
SURREY	793702	5742	784542	6032	1.17
HERTFORDSHIRE	723683	3622	677935	3670	6.75
BUCKINGHAMSHIRE	700343	2778	648867	2882	7.93
BRIGHTON AND HOVE	588432	516	588451	544	-0.00
OXFORDSHIRE	562062	2584	532447	2760	5.56
WOKINGHAM	554932	1009	521463	1016	6.42
READING	539202	246	472896	280	14.02
BRACKNELL FOREST	531927	482	466597	495	14.00
WEST BERKSHIRE	527765	782	503797	831	4.76
WEST SUSSEX	491486	4233	468644	4466	4.87
BATH AND NORTH EAST SOMERSET	465811	614	490501	625	-5.03
HAMPSHIRE	463830	7483	440626	7603	5.27
POOLE	459666	1213	426080	1172	7.88

Alla on kuva Windsorin ja kaikkien alueiden hintajakaumista. Tarkastelu on tässä rajattu alle kahden miljoonan punnan vanhoihin omakotitaloihin. Jakaumat ovat positiivisesti vinoja, joten ei ole yllätys, että kummankin jakauman keskiarvo ylittää niiden mediaaniarvon.

Alla on vielä karttakuva, josta selviävät asuntojen keskihinnat alueittain. Tätä kuvaa varten on ladattava Englannin koordinaatit ja yhdistettävä ne alkuperäiseen dataan. Tämän jälkeen kuva voidaan muodostaa esim. qplotilla tai ggplotilla. Alla olevan kuvan tein qplotilla.

R toimii suurilla aineistolla suhteellisen hyvin, mutta selkeä yläraja alkaa tulla vastaan, kun datan koko ylittää 3 gigaa. Ongelmia suurten aineistojen kanssa tulee varsinkin, jos tavoitteena on mallintaa dataa esim. regressiomalleilla. Tällöin laskemisajat voivat olla jopa tunteja. R:lle on kehitetty pakkaus nimeltä SparkR, joka pyrkii parantamaan R:n toimintaa Big datan kanssa. Pakkaus yhdistää Sparkin ja R:n, eli Sparkia käytetään R:n kautta. Kirjoittelen SparkR:sta jatkossa lisää.

Ristiintaulukointi ja CLT-simulointi

Tässä tekstissä käsittelen ristiintaulukointia sekä muita hyödyllisiä datan tarkastelumenetelmiä.
Ristiintaulukointi saadaan tehtyä table-komennolla. Komento laskee kahdelle luokitellulle muuttujalle niiden väliset frekvenssit. Aloitetaan simuloimalla yritysdataa, jossa on mukana toimialaluokitus, tieto mahdollisesta konkurssista sekä nettotulosprosentti.

R-koodi:

toimiala <- sample(letters[1:4],1000,replace=TRUE)
konkurssi <- rbinom(1000,1,0.02)
nettotulos <- rnorm(500,mean=0,sd=5)+runif(500,-5,25)
data <- data.frame(toimiala,konkurssi,nettotulo

Aloitetaan datan tutkiminen toimialaluokitus-muuttujasta katsomalla eri toimialaluokkien osuudet.

table(data$toimiala)

toimiala
  a   b   c   d
241 267 244 248

Ratkaistaan myös prosenttiosuudet

prop.table(table(data$toimiala))

toimiala
    a     b     c     d
0.241 0.267 0.244 0.248

Visualisoidaan barplotilla ja piechartilla:

aa <-barplot(prop.table(table(toimiala)),col=c("5","6","2","3"),ylim=c(0,1),main="Eri toimialaluokkien osuudet kaikista")
text(aa,0.15,round(prop.table(table(toimiala)),2),cex=1.3,pos=3)

library(plotrix)
tabl <- prop.table(table(toimiala))
lbl <- paste(names(tabl), "\n", tabl, sep="")
pie3D(tabl,labels=lbl,explode=0.05,main="Eri toimialojen osuudet kaikista")

Eri toimialoista on siis lähes saman verran havaintoja.

Katsotaan seuraavaksi barplot konkurssi-muuttujasta.

a <- barplot(prop.table(table(konkurssi)),col=c("10","1"),ylim=c(0,1),main="Konkurssit")
text(a,0.2,round(prop.table(table(konkurssi)),3),cex=1.1,pos=3)
legend("topright",legend=c("Konkurssi","Ei konkurssia"),pch=c(19,19),col=c(10,1))

Konkurssiin on mennyt datan tuhannesta yrityksestä 2.5 prosenttia. Aloitetaan nyt ristiintaulukointi ja tarkastellaan konkursseja toimialoittain. Ensin table-komennolla:

table(konkurssi,toimiala)

         toimiala
konkurssi   a   b   c   d
        0 233 258 240 244
        1   8   9   4   4

prop.table(table(konkurssi,toimiala),margin=2)

        toimiala
konkurssi          a          b          c
        0 0.96680498 0.96629213 0.98360656
        1 0.03319502 0.03370787 0.01639344
         toimiala
konkurssi          d
        0 0.98387097
        1 0.01612903

Tehdään vielä barplot kuva konkurssien osuuksista eri toimialoilla

a <- barplot(prop.table(table(konkurssi,toimiala),2),beside=T,yaxt="n",ylab="%",xlab="Toimiala") # Tehdään pylväskuvio
legend("topright",legend=c("Ei konkurssia","Konkurssi"),bg="white",pch=c(19,19),col=c(1,"grey"))
axis(2,at=seq(0,1,by=0.1))
text(a, 0.1, round(prop.table(table(konkurssi,toimiala)),2),cex=1,pos=3)

Testataan prop.testilla, onko eri toimialojen välillä tilastollisesti merkittävää eroa konkurssien osuuksissa. Menetelmällä testataan nollahypoteesia, jonka mukaan osuudet ovat kaikissa toimialoissa samat. Nollahypoteesia ei kuitenkaan hylätty.

• H₀: The proportion of cases is the same in each age group: p₁ = p₂ = p₃ = p₄ = p₅ = p₆
• H1: The proportion of cases is not the same in each age group: at least one p_i is different from the others

yhteensävektori = tapply(data$konkurssi, data$toimiala, length)
konkurssejavektori = tapply(data$konkurssi[data$konkurssi==1],data$toimiala[data$konkurssi==1], length)
prop.test(konkurssejavektori,yhteensävektori)


Katsotaan seuraavaksi keskimääräiset nettotulokset eri toimialoilla ja eri konkurssi-muuttujan arvoilla.

tapply(data$nettotulos,list(data$konkurssi,data$toimiala),mean)

      a        b         c         d
0 9.963086 9.310432 10.066982  9.532678
1 4.142009 5.814025  1.901524 18.058030

Toimialoilla a,b ja c konkurssiin menneillä yrityksillä oli alempi nettulosprosentti.

z <- tapply(data$nettotulos,list(data$konkurssi,data$toimiala),mean)
a=barplot(z,beside=T,main="Nettotulosprosentti toimialan ja konkurssi-statuksen mukaan",col=c(12,14),ylim=c(0,20),args.legend=list(bg="white",x=12,y=21),legend.text=c("Ei konkurssia","Konkurssi"))
text(a,0.4, round(z,1),cex=1.3,pos=3)

Central Limit Theorem -simulointi

Lopuksi näytän vielä, kuinka Central Limit Theoremin takia voimme yleistää otoksesta saadut keskiarvot koko populaatioon riippumatta populaation jakaumasta. Saamamme luottamusvälit ovat siis valideja vaikka esim. nettotulosprosentin jakauma olisi vinoutunut. Testataan tätä R:n For-loopin avulla. Muodostetaan ensin Poisson-jakautunutta dataa, joka edustaa populaation jakaumaa ja on tarkoituksella kaukana normaalijakaumasta. Tästä populaatiosta otamme sitten otoksia ja testaamme, osuuko populaation keskiarvo otoksesta saadun keskiarvon luottamusvälille. Lasketaan otoksen keskiarvoille luottamusväli 95%:n mukaan, eli simuloinnista voimme odottaa, että n.5%  kerroista populaation keskiarvo ei osu luottamusvälillemme. Alla kuva jakaumasta.

R-koodi:

a <- rpois(100000,1)
hist(a, main="Poisson-jakauma",col="5",xlab="")
a <- rpois(100000,1)


hist(a)
a <- data.frame(a)
summary(a)


tulosvektori <- vector("numeric",10000)  #Vektori,johon loopin tulokset tulevat


se <- function(x){                  #Standard errorin funktio
sd(x)/sqrt(length(x))
}



for (i in 1:10000){                                      
set.seed(i)
d <- sample(nrow(a),0.02*nrow(a))
b <- a[d,]
tulosvektori[i] <- ifelse(mean(a[,1])>=mean(b)-1.96*se(b)& mean(a[,1]) <=mean(b)+1.96*se(b),1,0)
}

prop.table(table(tulosvektori))

tulosvektori
    0     1
0.047 0.953

Populaation oikea keskiarvo osui 10000:sta lasketusta luottamusvälistä 95.3%:iin, eli 95%:n luottamusväli piti hyvin paikkansa.

Vastaava tulos saatiin myös alla olevalla beta-jakaumalle.

Nämä tulokset ja luottamusvälin toiminta perustuvat siihen, että otoksista lasketut keskiarvot noudattavat normaalijakaumaa. Tämä voidaan nähdä katsomalla otoskeskiarvojen jakaumaa.

R-koodi:

for (i in 1:10000){
set.seed(i)
d <- sample(nrow(a),0.02*nrow(a))
b <- a[d,]
vektori[i] <- mean(b)
}