####################
# R-harjoitukset 2 #
# viikko 37, 2007  #
####################

############
# teht. 1. #
############

# Aineiston lukeminen ja nimeäminen:

# useita vaihtoehtoja: 
# onn <- read.table("C:\\ . . .\\kaivos_onn.dat",header=T,skip=5)
# onn <- read.table("C:/ . . . /kaivos_onn.dat",header=T,skip=5)
# Jos hakemistopolku on määritelty jo aiemmin, ei sitä enää tarvita
# onn <- read.table("kaivos_onn.dat",header=T,skip=5)
# Polku voi viitata myös www-osoitteeseen:

onn <- read.table("http://mtl.uta.fi/tilasto/mtt06/Datat/kaivos_onn.dat",header=T,skip=5)
# huomaa taas lisämääreet: header=T ja skip=5

# Muuttujiin voidaan viitata $-merkin avulla, esim. tässä
onn$days

# Liitetään muuttujat työskentelytilaan, jonka jälkeen aineiston
# muuttujiin voidaan viitata suoraan nimillä

attach(onn)
names(onn)

# a)

# Lasketaan aineistosta perustunnuslukuja:
mean(days)
sd(days)
median(days)
min(days); max(days)
summary(days)

# piirretään kuvaaja:
hist(days)


# b)

# (i) #

# histogrammiin omat luokkarajat ja piirretään uusi kuvaaja:
br <- c(seq(-0.5,399.5,50),seq(499.5,999.5,100),1999.5)
hist(days,br)

# vastaava frekvenssitaulukko:
frekvenssit <- table(cut(days,br)) ; frekvenssit

# (ii) #

# Piirretään useampi kuvaaja samaan kuvaikkunaan, jolloin vertailu helpottuu:
par(mfrow=c(2,2))
hist(days[1:56],br) # histogrammi 56 ensimmäisestä havainnosta, todelliset luokat
hist(days[57:109],br) # histogrammi 53 viimeisestä havainnosta, todelliset luokat
hist(days,br) # vertailun vuoksi koko aineisto, todelliset luokat

# (iii) #

# huomataan, että kaivosten turvallisuus on parantunut, koska jälkimmäisellä jaksolla
# suhteellisesti enemmän pitkiä aikoja onnettomuuksien välillä

# vertaa "huonoon" tapaan:
par(mfrow=c(2,2))
hist(days[1:56])
hist(days[57:109])
hist(days)

# palautetaan ns. tavallinen asetus kuvaikkunaan (yksi kuvaaja ikkunassa):
par(mfrow=c(1,1))


# c)

# empiirinen kertymäfunktio:

Fn <- ecdf(days)
plot(Fn, verticals= TRUE, do.p = FALSE)
# tai
plot.ecdf(days, verticals= TRUE, do.p = FALSE)


# d)

# (i) #

# empiirisen kertymäfunktion arvon määrääminen pisteessä 60
Fn(60)
# tai
length(which(days<=60))/length(days)

# (ii) #

# P(149.5,199.5) ekf:n avulla
Fn(199.5)-Fn(149.5)
# tai
(length(which(days<=199.5))-length(which(days<=149.5)))/length(days)


############
# teht. 2. #
############
# normaalijakauma

# a)

par(mfrow=c(2,1)) # piirretään kaksi seuraavaa kuviota samaan ikkunaan
hist(rnorm(100)) # R olettaa N(0,1)-jakauman automaattisesti
hist(rnorm(100,50,5)) # sama asia on antaa komento: hist(rnorm(100,mean=50,sd=5))

# b)

dnorm(0) # N(0,1)-jakauman tf-arvo, kun X=0
dnorm(20,30,8) # N(30,8^2)-jakauman tf-arvo, kun X=20
dnorm(20,mean=30,sd=8) # antaa saman kuin edellä

par(mfrow=c(1,2))
plot(seq(-4,4,0.1),dnorm(seq(-4,4,0.1)),type="l") # N(0,1)-jakauman tapaus
plot(seq(0,60,0.1),dnorm(seq(0,60,0.1),30,8),type="l") # N(30,8^2)-jakauman tapaus

# ekstra: lisätään viimeksi tehtyyn kuvaan toinen normaalijakauman tiheysfunktio, 
# jossa odotusarvo on sama ja keskihajonta on 12. Kuvaaja katkoviivoin:
curve(dnorm(x,30,12),0,60,add=T,lty=2)

# c)

pnorm(0) # N(0,1)-jakauman kf-arvo, kun X=0
pnorm(20,30,8) # N(30,8^2)-jakauman kf-arvo, kun X=20

par(mfrow=c(1,2))
plot(seq(-4,4,0.1),pnorm(seq(-4,4,0.1)),type="l")
plot(seq(0,60,0.1),pnorm(seq(0,60,0.1),30,8),type="l")

# ekstra: lisätään viimeksi tehtyyn kuvaan toinen normaalijakauman kertymäfunktio, 
# jossa odotusarvo on sama ja keskihajonta on 12. Kuvaaja katkoviivoin:
lines(seq(0,60,0.1),pnorm(seq(0,60,0.1),30,12),lty=2)


# d)

qnorm(0.95) # N(0,1)-jakauman tapauksessa 95 %:n kvantiili
qnorm(0.95,3,1) # N(3,1)-jakauman tapauksessa 95 %:n kvantiili


############
# teht. 3. #
############
# binomijakauma

# a)

sum(dbinom(46:54,100,0.5))
# tai
pbinom(54,100,0.5)-pbinom(45,100,0.5)

# b)

rbinom(10,50,0.2)


############
# teht. 4. #
############
# kirjan luvun 1 harjoitus 4.
# harhattoman nopan heitto, sample-funktio
?sample

br <- seq(0.5,6.5,1) # mielekäs luokkajako

# esimerkki aluksi: heitetään 60 kertaa harhatonta noppaa ja tutkitaan sitä:
y <- sample(1:6,60,replace=TRUE,rep(1/6,6)) 
table(cut(y,br))/60 # suhteelliset frekvenssit taulukossa
par(mfrow=c(1,2)) # saadaan samaan kuvioon kaksi histogrammia
hist(y, br, freq = TRUE) # frekvenssihistogrammi
hist(y, br, freq = FALSE) # suhteellisten frekvenssien histogrammi

# tehdään heittoja 60, 120, 240, 480, 960 ja 2000:
y60 <- sample(1:6, 60, replace=TRUE, rep(1/6,6))
y120 <- sample(1:6, 120, replace=TRUE, rep(1/6,6))
y240 <- sample(1:6, 240, replace=TRUE, rep(1/6,6))
y480 <- sample(1:6, 480, replace=TRUE, rep(1/6,6))
y960 <- sample(1:6, 960, replace=TRUE, rep(1/6,6))
y2000 <- sample(1:6, 2000, replace=TRUE, rep(1/6,6))

# eri silmälukujen suhteelliset frekvenssit taulukoissa eri heittosarjoissa:
table(cut(y60,br))/60 
table(cut(y120,br))/120
table(cut(y240,br))/240
table(cut(y480,br))/480
table(cut(y960,br))/960
table(cut(y2000,br))/2000

# suhteellisten frekvenssien histogrammat:
par(mfrow=c(3,2))
hist(y60, br, freq = FALSE)
hist(y120, br, freq = FALSE)
hist(y240, br, freq = FALSE)
hist(y480, br, freq = FALSE)
hist(y960, br, freq = FALSE)
hist(y2000, br, freq = FALSE)

# Suhteelliset frekvenssit lähenevät lukua 1/6, kun n kasvaa.

# palautetaan kuvaajaikkuna tavalliseksi:
par(mfrow=c(1,1))

# jos haluaisi tehdä pylväsdiagrammin suhteellisista frekvensseistä histogrammin sijaan,
# se onnistuu barplot-funktion avulla, esim. 60 heiton tapauksessa:
barplot(table(cut(y60,br))/60, names=c("1","2","3","4","5","6"))


############
# teht. 5. #
############
# kirjan luvun 1 harjoitus 9.
# kahden harhattoman nopan heitto samanaikaisesti

# b)

# 100 heittoa:
n <- 100

# todennäköisyydet:
p <- c(1,2,3,4,5,6,5,4,3,2,1)/36

# tehdään heitot ja kerätään tulokset y-muuttujaan:
y <- sample(2:12,n,replace=TRUE,p)

br <- seq(1.5,12.5,1) # mielekäs luokkajako

# Havaitut frekvenssit:
hist(y,br) 
table(cut(y,br))

# Odotetut frekvenssit (n*p):
100*p


############
# teht. 6. #
############
# kolmen nopan heitto samanaikaisesti

tulos <- expand.grid(noppa1=1:6, noppa2=1:6, noppa3=1:6)
# saadaan matriisi, jonka kullakin rivillä on yksi tulosvaihtoehto

dim(tulos) # funktiolla dim saadaan matriisin rivien ja sarakkeiden lukumäärät
attach(tulos) # liitetään muuttujien nimet työtilaan (noppa1,noppa2,noppa3)

# kaikki sellaiset rivit, joissa kaikki silmäluvut samoja (suotuisat tapaukset):
n<-tulos[noppa1==noppa2&noppa2==noppa3,]
n

dim(n)[1] # suotuisten tapausten lkm (n-matriisin ensimmäinen dimensio eli rivit)
dim(tulos)[1] # kaikkien tapausten lkm (tulos-matriisin ensimmäinen dimensio eli rivit)

# kysytty todennäkäisyys:
dim(n)[1]/dim(tulos)[1]

# kysytty todennäköisyys osamääränä (tarkka arvo)
library(MASS)
fractions(dim(n)[1]/dim(tulos)[1])