### Eksempler Kap. 9, spesielt tester på data.

## 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)
dim(tips)
names(tips)

# Lager et sammendrag av dataene
summary(tips) 
boxplot(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.tips = length(tips$TIPpc)
tips.mean = mean(tips$TIPpc)
tips.sd = sd(tips$TIPpc)
z.obs = (tips.mean-15)/(tips.sd/sqrt(n.tips))

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

# Men vi bør vel gjøre t-korreksjon for estimert sigma, t.obs blir lik z.obs, 
# men vi forkaster først når t.obs>t_(0.05)
t_0.05=qt(0.95,n.tips-1);t_0.05
# altså t_0.05=1.67 < t.obs=4.21, forkaster fortsatt

# t-testen (med konfidensintervaller) gjøres oftest ved kommando t.test.
t.test(tips,mu=15)


## Styrkeberegninger med t-test
## Eksempel 9.14 (Voltage drop)

# 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")

# Figur: Styrke tosidige t-tester med ulike n
x=seq(-4,4,0.0001)
par(mfrow=c(1,1))
plot(x,power.t.test(n=100,delta=x,sd=10,type="one.sample",alternative="two.sided")$power,type="l",ylab="Power",xlab="delta")
lines(x,power.t.test(n=50,delta=x,sd=10,type="one.sample",alternative="two.sided")$power,col=2)
lines(x,power.t.test(n=20,delta=x,sd=10,type="one.sample",alternative="two.sided")$power,col=3)
legend(-3,0.95,col=c(1,2,3),c("n=100","n=50","n=20"),bty="n",lty=rep(1,3))

## Eksempel 9.12 Forurensning langs veien

# 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

#t-tetthet med 8 frihetsgrader
x=seq(-4,4,0.01)
plot(x,dt(x,8),type="l",ylab="t-tetthet, df=8")
abline(v=qt(0.995,n-1),col=2)
abline(v=-qt(0.995,n-1),col=2)
abline(v=t.obs)

# P-verdi for svevstøvdataene, P(ny t.obs>3.59, H0: mu=44, n=9)
2*pt(-t.obs,df=n-1)

# Framover vil vi bruke kommandoen t.test. Får da t.obs, p-verdi samt konfidensintervall direkte
t.test(d50,mu=44,alternative="two.sided")

# Tilsvarende t-test for tipsdataene (Eks. 9.13)
t.test(tips,mu=15,alternative="greater")

# (Får da i tillegg til t.obs og p-verdi også
# et "ensidig konfidensintervall", se s. 460 i læreboka)

# Ved å teste mot nedre grense i dette intervallet får vi
t.test(tips,mu=16.803,alternative="greater")
# altså p-verdi lik nivået 5%
# og mot gjennomsnittet
t.test(tips,mu=17.98614,alternative="greater")
# altså p-verdi=50%


## Forsoek med salve
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.05
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%

# P-verdi naar vi Y = 60
1-pnorm((p.hat-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(Y,n,p=p.0,alternative="greater",correct=F)
# Merk at R gir oss z.obs^2 = X-squared

# Vi kan ogsaa utfoere testen med korreksjon for kontinuitet (jf s. 223-225 i laereboka)
prop.test(Y,n,p=p.0,alternative="greater")


# Eksakt test
1-pbinom(Y-1,n,p.0) #
# Med korreksjon for kontinuitet faar vi nesten samme P-verdi som med den eksakte testen

# Den eksakte testen kan alternativ gjøres med binom.test (som også gir ensidig KI, her)
binom.test(Y,n,alternative="greater")

# Eventuelt for alternativ ulik p0=0.5
binom.test(Y,n,alternative="two.sided")
# som også gir eksakt konfidensintervall verdt å sammenligne med prop.test KI
prop.test(Y,n,0.5,alternative="two.sided")
# som altså er ganske like,




# Styrkefunksjon for mu med tosidig alternativ mu ulik mu0=2 kjent sigma=1, nivå 0.05
z0.05=qnorm(0.975);mu0=2;sigma=2
styrke=function(m,n) 1-pnorm(z0.05+(mu0-m)*sqrt(n)/sigma)+pnorm(-z0.05+(mu0-m)*sqrt(n)/sigma)
m=seq(1,3,0.001)
plot(m,styrke(m,10),type="l",ylab="styrke",ylim=c(0,1))
for (i in 2:8) lines(m,styrke(m,10*i),col=i)
abline(h=0.05,lty=2)
#abline(h=0.8,lty=3)
legend(1.5,1.05,c("n=10","n=20","n=30","n=40","n=50","n=60","n=70","n=80"),lty=rep(1,8),col=1:8,bty="n")


# Dobbel eksponensiell fordeling f(x;theta)=exp(-|x-theta|) 
theta=2
x=theta+seq(-4,4,0.001)
plot(x,0.5*exp(-abs(x-theta)),type="l",ylab="tetthet",main="Dobbel eksponensiell fordeling")
text(-1,0.95,"theta=2")  

# Simulerer n=3 verdier fra dobbel eksponensiell fordeling og plotter -loglik
n=3
x=theta+rexp(n)*(2*rbinom(n,1,0.5)-1)
loglik=numeric(0)
thetakand=seq(min(x)-1,max(x)+1,0.001)
for (i in 1:length(thetakand)) loglik[i]=sum(abs(x-thetakand[i]))
plot(thetakand,loglik,type="l",xlab="theta",main="-log-lik dobbel eksp")
sort(x)
for (i in 1:n) abline(v=x[i],lty=3)

# Vil sjekke -2ln(LR) tiln kji-kvadrat, df=1
n=100
theta=theta0=1
minus2lnLR=numeric(0)
for (sim in 1:1000){
x=theta+rexp(n)*(2*rbinom(n,1,0.5)-1)
thetahat=median(x)
minus2lnLR[sim]=2*(sum(abs(x-theta0))-sum(abs(x-thetahat)))
}
hist(minus2lnLR,prob=T,nclass=50,xlim=c(0,8))
v=seq(0.001,15,0.001)
lines(v,dchisq(v,1))
legend(1.5,1,"kji-kvadrat-tetthet, df=1",lty=1,bty="n")

## Test for om nedboersdata 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]
boxplot(regn)
hist(regn)

# Beregner momentestimatene (DBC, s. 417)

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))
# Testobservator og forkastningsgrense
show(c(minus.2lnLR,qchisq(0.95,1)))
# P-verdi
1-pchisq(minus.2lnLR,1)



## 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("mu = 30"),ylim=c(0,9))
hist(p.val.1,probability=T,xlab="P-verdi",main=paste("mu = 31"),ylim=c(0,9))
hist(p.val.2,probability=T,xlab="P-verdi",main=paste("mu = 32"),ylim=c(0,9))

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


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