# KOMMANDOER TIL FORELESNINGENE FOR UKE 16 

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

# EKSEMPEL: RUTHERFORD OG GEIGERS STRÅLINGSDATA
# (test om Poisson-fordeling)

# Leser inn dataene:
alfa.partikler=read.table("http://www.uio.no/studier/emner/matnat/math/STK2120/v13/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 når vi slår 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 på N-ene

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

# Finner ML-estimatoren ved å bruke kommandoen "optimize" (som finner 
# maksimum av en funksjon av én 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 nøyere en seg med å bestemme ML-estimatoren for de 
# opprinnelige dataene (X-ene); # jf sidene 721-722 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 722 i D&B 

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


# EKSEMPEL: FØDSELSVEKT FOR GUTTER (test av normalitet)

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

# Ser på de 338 guttene som har svangerskapslengde 40 uker
vekt=fvekt$vekt[(fvekt$varighet==40)&(fvekt$kjonn==1)]

# Histogram
hist(vekt,breaks=c(1500,2500,3000,3500,4000,4500,5500))

# Observerte antall i intervallene gitt ved "breaks":
N=hist(vekt,breaks=c(1500,2500,3000,3500,4000,4500,5500),plot=F)$counts

# Vil finne ML-estimator for my og sigma basert på N-ene
# Definerer første en funksjon som beregner minus log-likelihood:
negloglik=function(par){
my=par[1];s=par[2]
prob=pnorm(c(2500,3000,3500,4000,4500,Inf),my,s)-
     pnorm(c(-Inf,2500,3000,3500,4000,4500),my,s)
loglik=sum(N*log(prob))
-loglik}

# Finner ML-estimatoren ved å bruke kommandoen "optim" (som finner 
# maksimum av en funksjon av flere variabler):
par=optim(c(3200,400),negloglik)$par
my=par[1]
sigma=par[2]

# Forventete antall:
n=length(vekt)
prob=pnorm(c(2500,3000,3500,4000,4500,Inf),my,sigma)-
     pnorm(c(-Inf,2500,3000,3500,4000,4500),my,sigma)
E=n*prob

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

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

# ML-estimatoren for de opprinnelige dataene (X-ene) er gitt ved:
my.x=mean(vekt)
sigma.x=sd(vekt)*sqrt((n-1)/n)

# Det gir forventete antall hvis normalfordelt:
prob.x=pnorm(c(2500,3000,3500,4000,4500,Inf),my.x,sigma.x)-
     pnorm(c(-Inf,2500,3000,3500,4000,4500),my.x,sigma.x)
E.x=n*prob.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 3 og 5 frihetsgrader; jf. side 722 i D&B 


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

# EKSEMPEL: POLITISKE MENINGSMÅLINGER

# Leser inn data:
parti=c('RV','SV','Ap','Sp','MDG','V','KrF','H','FrP','Andre')
mars=c(15, 33, 201, 39, 8, 35, 33, 227, 130, 4)
april=c(9, 31, 218, 35, 6, 29, 33, 244, 113, 10)

# Oppslutning om partiene (prosent)
mars.prosent=round(100*mars/sum(mars),1)
april.prosent=round(100*april/sum(april),1)
rbind(parti,mars.prosent,april.prosent)

# Vi finner de forventete verdiene:
n1=sum(mars)
n2=sum(april)
n=n1+n2
E.mars=n1*(mars+april)/n
E.april=n2*(mars+april)/n

# Sammenligning av observerte og forventete verdier
rbind(mars,E.mars)
rbind(april,E.april)

# Beregning av kji-kvadrat observatoren med p-verdi
N=rbind(mars,april)
E=rbind(E.mars,E.april)
X2=sum((N-E)^2/E)
1-pchisq(X2,9)

# Vi kan også bruke kommandoen "chisq.test":
parti.test=chisq.test(rbind(mars,april))
parti.test

# Forventede verdier under nullhypotesen:
parti.test$expected