# KOMMANDOER TIL FORELESNINGENE FOR UKE 16 

# ========================================


# EKSEMPEL: KRYSSING AV TOMATER (eksempel 13.1 i D&B)
# En forklaring av dihybrid krysning er gitt her:
# https://en.wikipedia.org/wiki/Dihybrid_cross

# Leser inn dataene fra tabell 13.2:
N=c(926,288,293,104)

# Sannsynligheter under nullhypotesen om Mendelsk nedarving:
p0=c(9/16,3/16,3/16,1/16)

# Forventete antall:
n=sum(N)
E=n*p0

# Tabell over observerte og forventete antall:
rbind(N,E)

# Kji-kvadrat observator:
X2=sum((N-E)^2/E)

# P-verdi
1-pchisq(X2,3)

# Vi kan bruke kommandoen "chisq.test" til  gjre beregningene:
X2test=chisq.test(N,p=p0)
X2test
rbind(X2test$observed,X2test$expected)

# ================================


# EKSEMPEL: TID FOR START AV FOEDSELSVEER (eksempel 13.3 i D&B)

# Har tiden da veene starter for n=1186 foedsler.
# Spoersmaalet er om tidene fordeler seg uniformt over doegnet.

# Antall fdsler i timesintervaller etter midnatt
N=c(52,73,89,88,68,47,58,47,48,53,47,34,21,31,40,24,37,31,47,34,36,44,78,59)

# Sannsynligheter under nullhypotesen:
p0=rep(1/24,24)

# Bruke kommandoen "chisq.test" til  gjoere beregningene:
X2test=chisq.test(N,p=p0)
X2test
rbind(X2test$observed,X2test$expected)



 EKSEMPEL: FORDELING AV BARNAS KOEJNN I EN FAMILIE

# Data fra Helge Brunborg: Blir det gutt eller jente?
# Samfunnspeilet nr 3, 2009
# http://www.ssb.no/a/samfunnsspeilet/utg/200903/ssp.pdf

# Ser paa kjoennsfordelingen for de tre foerste barna for 368445
# kvinner foedt 1935-91 som har faatt minst tre barn (tabell 3).
# Vi teste nullhypotesen om at antall gutter er binomisk fordelt

# Leser inn antall kvinner med 0, 1, 2 og 3 gutter:
N=c(43941,130883,139951,53670)

# Antall gutter (G) og jenter (J):
g=0:3
G=sum(g*N)

# Andel gutter:
n=sum(N)
tetta=G/(3*n)

# Forventete antall hvis binomisk fordelt:
E=n*dbinom(g,3,tetta)

# Tabell over observerte og forventete antall:
rbind(N,E)

# Kji-kvadrat observator:
X2=sum((N-E)^2/E)

# Det er et signifikant avvik fra binomisk fordeling


# ========================================

# EKSEMPEL: RUTHERFORD OG GEIGERS STRAALINGSDATA
# (test om Poisson-fordeling)

# Leser inn dataene:
alfa.partikler=read.table("http://www.uio.no/studier/emner/matnat/math/STK2120/v16/alfa-partikler.txt", header=T)

# Antall intervaller
n=sum(alfa.partikler$obs)

# Stolpediagram av dataene:
barplot(alfa.partikler$obs,names.arg=0:15)

# Observerte antall nr vi slaar sammen intervallene 
# med 11 eller flere alfa-partikler:
N=c(alfa.partikler$obs[1:11],sum(alfa.partikler$obs[12:16]))

# Vil finne ML-estimator for lambda basert paa N-ene

# Definerer en funksjon som beregner log-likelihooden:
loglik=function(lam) 
sum(N[1:11]*log(dpois(0:10,lam)))+N[12]*log(1-ppois(10,lam))

# Finner ML-estimatoren ved  bruke kommandoen "optimize" (som finner 
# maksimum av en funksjon av en variabel over et intervall):
lambda.n=optimize(loglik,c(0,10),maximum=T)$maximum

# Beregner forventete antall hvis Poisson fordelt:
E=n*c(dpois(0:10,lambda.n),1-ppois(10,lambda.n))

# Tabell over observerte og forventete antall:
rbind(N,E)

# Kji-kvadrat observator og P-verdi
X2=sum((N-E)^2/E)
1-pchisq(X2,10)

# Vanligvis noeyere en seg med  bestemme ML-estimatoren for de 
# opprinnelige dataene (X-ene); jf sidene 737-738 i D&B. 
# Her er ML-estimatoren basert X-ene lik antall partikler per intervall:
lambda.x=sum(alfa.partikler$ant*alfa.partikler$obs)/n

# Det gir forventete antall hvis Poisson fordelt:
E.x=n*c(dpois(0:10,lambda.x),1-ppois(10,lambda.x))

# Kji-kvadrat observatoren blir da:
X2.x=sum((N-E.x)^2/E.x)

# Denne har en fordeling som ligger mellom kji-kvadrat fordelingene
# med 10 og 11 frihetsgrader; jf. side 737 i D&B