### Eksempler Kap. 9

## Forsoek med salve

# Type I-feil for forskjellige verdier av k
n = 100
p.0 = 0.5
k = 53:62
typeI.feil = 1-pbinom(k-1,n,p.0)
rbind(k,round(typeI.feil,3))
plot(k,typeI.feil)
abline(h=0.1,lty=2)

# Type II-feil for forskjellige verdier av p naar k = 59
n = 100
p = seq(0.55,0.8,0.05)
k = 59
typeII.feil = pbinom(k-1,n,p)
rbind(p,round(typeII.feil,5))
plot(p,typeII.feil)

# Styrkefunksjon for p mellom 0 og 1
par(mfrow=c(1,2))
n=100
p=seq(0,1,0.01)
k=59
styrke = 1-pbinom(k-1,n,p)
plot(p,styrke,type="l")

# Styrkefunksjon for p mellom 0.4 og 0.8
n = 100
p = seq(0.4,0.8,0.01)
k = 59
styrke.p = 1-pbinom(k-1,n,p)
plot(p,styrke.p,type="l")
abline(v=0.5,lty=2,lwd=2)

## Eksempel 9.3

# Styrkefunksjon for forskjellige verdier av mu
n = 10
sigma = 3
mu = seq(20,35,0.1)
c = 27.5
styrke.mu = 1-pnorm((c-mu)/(sigma/sqrt(n)))
plot(mu,styrke.mu,type="l")
abline(v=25,lty=2,lwd=2)

## Eksempel 9.13

# Leser dataene
data.fil = "https://www.uio.no/studier/emner/matnat/math/STK1110/data/exmp09-13.txt"
tips = read.table(data.fil,header=T)

# Lager et sammendrag av dataene
summary(tips)

# Lager histogram og normalfordelingsplot
par(mfrow=c(1,2))
hist(tips$TIPpc,xlab="Prosent tips",main="Histogram av prosent tips") 
qqnorm(tips$TIPpc)
qqline(tips$TIPpc)

# Tester H_0: mu <= 15 mot H_a: mu > 15
n = length(tips$TIPpc)
tips.mean = mean(tips$TIPpc)
tips.sd = sd(tips$TIPpc)
z.obs = (tips.mean-15)/(tips.sd/sqrt(n))

# En test med signifikansnivaa 5% forkaster H_0 hvis z.obs >= z_(0.05) = 1.645. 
z.obs

# Altsaa forkaster vi nullhypotesen

## Eksempel 9.12

# Leser dataene
d50 = c(82.9, 56.8, 66.5, 49.4, 105.4, 79.5, 82.5, 50.7, 43.0)

# Tester H_0: mu = 44 mot H_a: mu != 44
n = length(d50)
d50.mean = mean(d50)
d50.sd = sd(d50)
t.obs = (d50.mean-44)/(d50.sd/sqrt(n))

# En test med signifikansnivaa 1% forkaster H_0 hvis t.obs <= -qt(0.995,n-1)
# eller t.obs >= qt(0.995,n-1). 
# Vi ser paa verdien av test-observatoren og den kritiske verdien
t.obs
qt(0.995,n-1)

# Vi ser at t.obs >= qt(0.995,n-1), saa vi forkaster H_0

## Eksempel 9.14

# Vi ser paa testing av H_0: mu <= 2.5 mot H_a: mu > 2.5
# Vi finner styrken til testen naar mu = mu_1 = 2.6 
# (I komandoen nedenfor er delta = mu_1 - mu_0 = 2.6-2.5 = 0.10)
power.t.test(n=10,delta=0.1,sd=0.1,power=NULL,type="one.sample",alternative="one.sided")

power.t.test(n=NULL,delta=0.1,sd=0.1,power=0.8,type="one.sample",alternative="one.sided")


# Sannsynligheten for feil av type II er C)n minus teststyrken

# Plott av styrkefunksjonen
delta.vec = seq(-0.05,0.2,0.01)
styrke = power.t.test(n=10,delta=delta.vec,sd=0.1,power=NULL,type="one.sample",alternative="one.sided")$power
plot(delta.vec,styrke,type="l")
abline(v=0,lty=2,lwd=2)

# Vi finner n slik at testen faar sannsynlighet for type II-feil 5%, og
# dermed styrke 95% naar mu = mu_1 = 2.6
power.t.test(n=NULL,delta=0.1,sd=0.1,power=0.95,type="one.sample",alternative="one.sided")

## Forsoek med salve (forts.)

n = 100
p.0 = 0.5
Y = 60

# Da np_0 = 50 > 10 og n(1-p_0) = 50 > 10, kan vi gjenta testen basert paa normalfordeling
alpha = 0.0225
p.hat = Y/n
z.obs = (p.hat-p.0)/sqrt(p.0*(1-p.0)/n)
z.alpha = qnorm(1-alpha)
c(z.obs,z.alpha)

# Da z.obs >= z.alpha, forkaster vi H_0 ved signifikansnivaa 5%

# Styrkefunksjon for p mellom 0 og 1
p=seq(0,1,0.01)
styrke.n = 1-pnorm((p.0-p+z.alpha*sqrt(p.0*(1-p.0)/n))/sqrt(p*(1-p)/n))
plot(p,styrke.n,type="l",ylab="styrke")

# Legger til styrkefunksjon fra eksakt test med signifikansnivaa 4.4%
k=59
styrke.b = 1-pbinom(k-1,n,p)
lines(p,styrke.b,col=2)

# Sannsynlighet for type II-feil naar p = 0.6
p.1 = 0.6
pnorm((p.0-p.1+z.alpha*sqrt(p.0*(1-p.0)/n))/sqrt(p.1*(1-p.1)/n))

# Finne n for at denne skal bli hoeyst beta = 0.1
beta = 0.1
z.beta = qnorm(1-beta)
n = ((z.alpha*sqrt(p.0*(1-p.0))+z.beta*sqrt(p.1*(1-p.1)))/(p.1-p.0))^2
n
  
# Vi runder opp til n = 211.

## Forsoek med salve (forts.)

# P-verdi naar vi x.obs = 60
n = 100
p.0 = 0.5
x.obs.1 = 60
1-pbinom(x.obs.1-1,n,p.0) # Eksakt test
1-pnorm((x.obs.1/n-p.0)/sqrt(p.0*(1-p.0)/n)) # Med normalfordelingstilnaerming

# P-verdi naar vi x.obs = 65
x.obs.2 = 65
1-pbinom(x.obs.2-1,n,p.0) # Eksakt test
1-pnorm((x.obs.2/n-p.0)/sqrt(p.0*(1-p.0)/n)) # Med normalfordelingstilnaerming

# Vi kan ogsaa bruke kommandoen "prop.test" til aa utfoere testen
help(prop.test)
prop.test(x.obs.1,n,p=p.0,alternative="greater",correct=F)
prop.test(x.obs.2,n,p=p.0,alternative="greater",correct=F)

# Vi kan ogsaa utfoere testen med korreksjon for kontinuitet (jf s. 188-190 i laereboka)
prop.test(x.obs.1,n,p=p.0,alternative="greater")
prop.test(x.obs.2,n,p=p.0,alternative="greater")

# Med korreksjon for kontinuitet faar vi nesten samme P-verdi som med den eksakte testen

## Eksempel 9.13 (forts.)

data.fil = "https://www.uio.no/studier/emner/matnat/math/STK1110/data/exmp09-13.txt"
tips = read.table(data.fil,header=T)
n = length(tips$TIPpc)
tips.mean = mean(tips$TIPpc)
tips.sd = sd(tips$TIPpc)
z.obs = (tips.mean-15)/(tips.sd/sqrt(n))

# Vi finner P-verdien ved kommandoen
1-pnorm(z.obs)

## Eksempel 9.12 (forts.)

# Leser dataene
d50 = c(82.9, 56.8, 66.5, 49.4, 105.4, 79.5, 82.5, 50.7, 43.0)

# Tester H_0: mu = 44 mot H_a: mu != 44
n = length(d50)
d50.mean = mean(d50)
d50.sd = sd(d50)
t.obs = (d50.mean-44)/(d50.sd/sqrt(n))

# Her blir P-verdien
2*(1-pt(abs(t.obs),n-1))

# Vi kan bruke kommandoen "t.test" til aa utfoere testen
t.test(d50,alternative="two.sided",mu=44)

## Eksempel 9.22

N = 10000
n = 4
sigma = 2

# For mu = mu.0
mu.0 = 30
p.val.0 = rep(0,N)
for (i in 1:N)
{
  x = rnorm(n,mu.0,sigma)
  t.obs = (mean(x)-mu.0)/(sd(x)/sqrt(n))
  p.val.0[i] = 1-pt(t.obs,n-1)
}

# For mu = mu.1
mu.1 = 31
p.val.1 = rep(0,N)
for (i in 1:N)
{
  x = rnorm(n,mu.1,sigma)
  t.obs = (mean(x)-mu.0)/(sd(x)/sqrt(n))
  p.val.1[i] = 1-pt(t.obs,n-1)
}

# For mu = mu.2
mu.2 = 32
p.val.2 = rep(0,N)
for (i in 1:N)
{
  x = rnorm(n,mu.2,sigma)
  t.obs = (mean(x)-mu.0)/(sd(x)/sqrt(n))
  p.val.2[i] = 1-pt(t.obs,n-1)
}

# Histogrammer av P-verdiene
par(mfrow=c(1,3))
hist(p.val.0,probability=T,xlab="P-verdi",main=paste("Histogram av P-verdier, mu = 30"),ylim=c(0,9))
hist(p.val.1,probability=T,xlab="P-verdi",main=paste("Histogram av P-verdier, mu = 31"),ylim=c(0,9))
hist(p.val.2,probability=T,xlab="P-verdi",main=paste("Histogram av P-verdier, mu = 32"),ylim=c(0,9))

# Andel P-verdier <= 5%
mean(p.val.0 <= 0.05)
mean(p.val.1 <= 0.05)
mean(p.val.2 <= 0.05)

## Test for om nedboersdataene er eksponensialfordelt
## Tester da H_0: alpha = 1 mot Ha_: alpha != 1

# Leser dataene fra fil. De gir nedboer (i mm) hver av de 92 dagene i mai-juli 2016
nedbor.blindern = read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/nedbor-blindern.txt",header=T)

# Tar bare med de 54 dagene det regnet:
regn = nedbor.blindern$nedbor[nedbor.blindern$nedbor > 0]


# Beregner momentestimatene

n = length(regn)
mom.alpha = mean(regn)^2/(mean(regn^2)-mean(regn)^2)
mom.beta = (mean(regn^2)-mean(regn)^2)/mean(regn)
c(mom.alpha,mom.beta)

# Bestemmer maximum likelihood estimatene 
loglik_gamma = function(alpha,beta,x)
{
  -n*alpha*log(beta)-n*lgamma(alpha)+(alpha-1)*sum(log(x))-sum(x)/beta
}

negloglik_gamma = function(log.par,x)
{
  alpha = exp(log.par[1])
  beta = exp(log.par[2])
  -loglik_gamma(alpha,beta,x)
}

# Maksimerer log-likelihooden (ved aa minimere den negative log-likelihooden numerisk):
ml.est = optim(log(c(mom.alpha,mom.beta)),negloglik_gamma,x=regn)  
ml.alpha = exp(ml.est$par[1])
ml.beta = exp(ml.est$par[2])
c(ml.alpha,ml.beta)  

# Merk at hvis alpha = 1, har vi eksponensialfordelingen, og 
# ML-estimatet av beta er da lik gjennomsnittet
ml.beta0 = mean(regn)
ml.beta0

# Beregner -2*ln(LR) og P-verdien for likelihood ratio testen
minus.2lnLR = 2*(loglik_gamma(ml.alpha,ml.beta,regn)-
                   loglik_gamma(1,ml.beta0,regn))
show(c(minus.2lnLR,qchisq(0.95,1)))
1-pchisq(minus.2lnLR,1)

# P-verdien < 5%, og vi forkaster altsC% H_0