# KOMMANDOER TIL FORELESNINGENE FOR UKE 16 

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


# 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 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 737 i D&B 

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


# EKSEMPEL: FoEDSELSVEKT 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 bare paa guttene 
vekt=fvekt$vekt[fvekt$kjonn==1]

# Histogram
intervaller=c(1000,2250,2500,2750,3000,3250,3500,3750,4000,4250,4500,4750,5500)
hist(vekt,breaks=intervaller)

# Observerte antall i de ulike intervallene:
N=hist(vekt,breaks=intervaller,plot=F)$counts

# Vil finne ML-estimator for my og sigma basert paa N-ene
# Definerer foerst en funksjon som beregner minus log-likelihood:
negloglik=function(par){
my=par[1];s=par[2]
prob=pnorm(c(2250,2500,2750,3000,3250,3500,3750,4000,4250,4500,4750,Inf),my,s)-
pnorm(c(-Inf,2250,2500,2750,3000,3250,3500,3750,4000,4250,4500,4750),my,s)
loglik=sum(N*log(prob))
-loglik}

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

# Forventete antall:
n=length(vekt)
prob=pnorm(c(2250,2500,2750,3000,3250,3500,3750,4000,4250,4500,4750,Inf),my.n,sigma.n)-
pnorm(c(-Inf,2250,2500,2750,3000,3250,3500,3750,4000,4250,4500,4750),my.n,sigma.n)
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,9)

# 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(2250,2500,2750,3000,3250,3500,3750,4000,4250,4500,4750,Inf),my.x,sigma.x)-
pnorm(c(-Inf,2250,2500,2750,3000,3250,3500,3750,4000,4250,4500,4750),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 9 og 11 frihetsgrader; jf. side 737 i D&B 



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

# EKSEMPEL: POLITISKE MENINGSMLINGER

# Mars: http://www.nrk.no/norge/mdg-gar-kraftig-tilbake-i-ny-maling-1.12844492
# April: http://www.nrk.no/norge/aprilmaling_-arsbeste-for-team-erna-1.12897320

# Leser inn data:
parti=c('Rodt','SV','Ap','Sp','MDG','KrF','V','H','FrP','Andre')
mars=c(15,39,306,64,28,46,43,232,180,15)
april=c(15,34,315,58,31,50,45,241,166,15)

# 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


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

# EKSEMPEL: KJOENN OG FORETRUKKET HAAND
# Vil teste nullhypotesen at kjoenn og foretrukket haand er uavhengige

# Leser inn data:
haand=c('Hoyre','Venstre','Begge')
menn=c(934,113,20)
kvinner=c(1070,92,8)

# Fordeling for hvert kjoenn (prosent)
menn.prosent=round(100*menn/sum(menn),1)
kvinner.prosent=round(100*kvinner/sum(kvinner),1)
rbind(haand,menn.prosent,kvinner.prosent)

# Tabell med observerte verdier
N=rbind(menn,kvinner)

# Vi finner de forventete verdiene:
n=sum(N)
N.rad=apply(N,1,sum)
N.kolonne=apply(N,2,sum)
E=rbind(rep(NA,3),rep(NA,3))
for (i in 1:2)
for (j in 1:3)
E[i,j]=(N.rad[i]*N.kolonne[j])/n

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

# Vi kan ogs bruke kommandoen "chisq.test" for  utfre testen:
haand.test=chisq.test(N)
haand.test

# Forventede verdier under nullhypotesen:
haand.test$expected