################################################
#### oppgave fra kap 13: nr 16,18,19,23,24: 
################################################

#### oppgave 16
N = c(6,24,42,59, 62, 44, 41, 14, 6, 2)
k = length(N)
n = sum(N)
theta.hat = sum(N*c(0:(k-1)))/n       #gjennomsnittet til alle observasjonene
m=length(theta.hat)

p.hat = function(i) { theta.hat^i*exp(- theta.hat)/factorial(i)}
E <- c()
for(i in 1:k) {E[i] = n*p.hat(i)   }

#combining de siste 2 gruppene til en gruppe. 
N = c(6,24,42,59, 62, 44, 41, 14, sum(6, 2))
k = length(N)
E <- c(E[1:8], sum(E[9:10]) )

show(cbind(N,E))
chitest = sum((N-E)^2/E)
chiAlpha = qchisq(0.05,k-1-m, lower.tail=FALSE)
Pvalue = 1-pchisq(chitest,k-1-m)
paste("Chisq=",chitest, ", ChiAlpha=", round(chiAlpha, 3), ", P-value=",round(Pvalue,3))



#### oppgave 18 (test av normalitet, s723 i D&B)
N = c(12,20,23,15,13)

# Bruker estimatene:
mean_hat = 0.173
sd_hat   = 0.066

#plot over normalfordelingen vi ?nsker ? test om passer til v?re observasjoner
plot(density(rnorm(10000, mean=mean_hat, sd=sd_hat)))
abline(v=c(0.1, 0.15, 0.2, 0.25), col="red")		#v?re intervaller/grupper

# Forventete antall:
prob	=pnorm(c(0.1, 0.15, 0.2, 0.25, Inf), mean=mean_hat, sd=sd_hat) - pnorm(c(-Inf,0.1, 0.15, 0.2, 0.25 ), mean=mean_hat, sd=sd_hat)
n	=sum(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,3)



#### oppgave 19
N = c(49,26,14,20,53,38)
k = length(N)
n = sum(N)

#finner theta.hat ved ? derivere likelihooden p? theta1 og theta2. 
#dette resulterer i 2 likninger med 2 ukjent som kan l?ses slik:  (matrise form a*theta.hat=b, der a er en 2*2 matrise og b en vektor med lengde 2) 
theta.hat <- solve(a=cbind(c(290,110), c(171, 229)), b=c(171, 110)) 

p.hat = c(theta.hat[1]^2,theta.hat[2]^2, (1 - theta.hat[1] - theta.hat[2])^2, 2*theta.hat[1]*theta.hat[2],2*theta.hat[1]*(1 - theta.hat[1] - theta.hat[2]), 2*theta.hat[2]*(1 - theta.hat[1] - theta.hat[2]) )
m=length(theta.hat)
E = n*p.hat
show(cbind(N,E))
chitest = sum((N-E)^2/E)
chiAlpha = qchisq(0.05,k-1-m, lower.tail=FALSE)
Pvalue = 1-pchisq(chitest,k-1-m)
paste("Chisq=",chitest, ", ChiAlpha=", round(chiAlpha, 3), ", P-value=",round(Pvalue,3))



##### oppgave 23
# Leser inn data:
N = rbind(c(73,49), c(71,64))
colnames(N) = c("95","94")
rownames(N) = c("vunnet", "tapt")

# Vi finner de forventete verdiene:
ni. <- rowSums(N)
n.j <- colSums(N)
E <- matrix(NA, ncol=dim(N)[2], nrow=dim(N)[1] )
for(i in 1:dim(N)[1]){
    for(j in 1:dim(N)[2]){
        E[i,j] <- n.j[j]*ni.[i]/sum(N)
    }
}

X2=sum((N-E)^2/E)
1-pchisq(X2,1)

# Vi kan ogs? bruke kommandoen "chisq.test":
parti.test=chisq.test(N)
parti.test

# Forventede verdier under nullhypotesen:
parti.test$expected





##### oppgave 24
N <- rbind(c(409,11,22,7,277), c(512,4,14,11,220))
ni. <- rowSums(N)
n.j <- colSums(N)
E <- matrix(NA, ncol=dim(N)[2], nrow=dim(N)[1] )
for(i in 1:dim(N)[1]){
    for(j in 1:dim(N)[2]){
        E[i,j] <- n.j[j]*ni.[i]/sum(N)
    }
}
chitest = sum((N-E)^2/E)
chiAlpha = qchisq(0.05,(dim(N)[2]-1)*(nrow=dim(N)[1]-1), lower.tail=FALSE)
Pvalue = 1-pchisq(chitest,(dim(N)[2]-1)*(nrow=dim(N)[1]-1))
paste("Chisq=",chitest, ", ChiAlpha=", round(chiAlpha, 3), ", P-value=",round(Pvalue,3))