### Eksempler Kap. 8

# Eks. 8.1 DBC Tensile styrke i titanlegeringer
n=31;xmean=1064;sigma=55
stderrmean=sigma/sqrt(n)
stderrmean
round(xmean+c(0,-1,1)*stderrmean,1)
round(xmean+c(0,-1,1)*2*stderrmean,1)
round(xmean+c(0,-1,1)*1.96*stderrmean,1)

# Dekningsplot (coverage)
plot(mean(rnorm(10,0,2))+1.96*c(-1,1)*2/sqrt(10),rep(1,2),type="l",ylim=c(0,100),xlim=c(-2.9,2.9),ylab="",xlab="")
for (i in 2:100) lines(mean(rnorm(10,0,2))+1.96*c(-1,1)*2/sqrt(10),rep(i,2))
abline(v=0)


## Konsentrasjonsmaalinger

# Leser inn konsentrasjonsmaalingene:
kons = c(2.52,2.49,2.62,2.45,2.56)
n = length(kons)
n

# Normalt qq-plott
qqnorm(kons)
qqline(kons)

# Beregner gjennomsnitt og empirisk standardavvik
gjsn = mean(kons)
st.av = sd(kons)
c(gjsn,st.av)

# Bestemmer oevre (5/2)%-kvantil i t-fordelingen:
t.kvant = qt(0.975,n-1)
t.kvant

# Beregner 95% konfidensintervall:
low = gjsn-t.kvant*st.av/sqrt(n)
up = gjsn+t.kvant*st.av/sqrt(n)
print(c(low,up),digits=4)

# Til sammenligning - kvantiler fra normalfordelingen, ignorerer usikkerhet i standardavvik
print(gjsn+c(-1,1)*1.96*st.av/sqrt(n),digits=4)

# Eks. 8.1 DBC Tensile styrke i titanlegeringer
n=31;xmean=1064;sigma=55
stderrmean=sigma/sqrt(n)
stderrmean
#t-korrigert intervall
round(xmean+c(0,-1,1)*qt(0.975,n-1)*stderrmean,2)
#(naivt) normkvantil intervall
round(xmean+c(0,-1,1)*1.96*stderrmean,2)
#t-persentil
qt(0.975,n-1)

# Dekningsgrad t-intervaller
simres=numeric(0)
n=5
for (sim in 1:100){
  x=rnorm(n,0,2)
  simres=rbind(simres,c(mean(x),sd(x)))
}
apply(simres,2,summary)

low=min(simres[,1]-qt(0.975,n-1)*simres[,2]/sqrt(n))
high=max(simres[,1]+qt(0.975,n-1)*simres[,2]/sqrt(n))
plot(simres[1,1]+c(-1,1)*qt(0.975,n-1)*simres[1,2]/sqrt(n),c(1,1),type="l",ylim=c(0,100),xlim=c(low,high),xlab="KI",ylab="simulering no.")
for (i in 2:100) lines(simres[i,1]+c(-1,1)*qt(0.975,n-1)*simres[i,2]/sqrt(n),c(i,i))
abline(v=0)
for (i in 1:100) lines(simres[i,1]+c(-1,1)*1.96*simres[i,2]/sqrt(n),c(i,i),col=2)

in_t_int=sum((simres[,1]-qt(0.975,n-1)*simres[,2]/sqrt(n)<0)*(simres[,1]+qt(0.975,n-1)*simres[,2]/sqrt(n)>0))
in_norm_int=sum((simres[,1]-1.96*simres[,2]/sqrt(n)<0)*(simres[,1]+1.96*simres[,2]/sqrt(n)>0))
c(in_t_int,in_norm_int)


# 95% prediksjonsintervall for en ny maaling:
low.pred=gjsn-t.kvant*st.av*sqrt(1+1/n)
up.pred=gjsn+t.kvant*st.av*sqrt(1+1/n)
print(c(low.pred,up.pred),digits=3)

# Bestemmer nedre oevre (5/2)%-kvantil i chi^2-fordelingen:
chi2.kvant.low = qchisq(0.025,df=n-1)
chi2.kvant.up = qchisq(0.975,df=n-1)
c(chi2.kvant.low,chi2.kvant.up)

print(c(st.av^2,st.av),digits=3)


# 95% konfidensintervall for variansen:
low.var=(n-1)*st.av^2/chi2.kvant.up
up.var=(n-1)*st.av^2/chi2.kvant.low
print(c(low.var,up.var),digits=3)

# 95% konfidensintervall for standardavviket:
print(sqrt(c(low.var,up.var)),digits=3)

# Simulerer 1000 normalfordelte datasett med n=10 og beregner gjennomsnitt og empirisk varians i hvert sett.
simres=numeric(0)
for (sim in 1:1000){
  x=rnorm(10,0,1)
  simres=rbind(simres,c(mean(x),var(x)))
}
plot(simres[,1],simres[,2],xlab="gjennomsnitt",ylab="empirisk varians",main="Normalfordelte data")
cor(simres[,1],simres[,2])
# Plott og korrelasjon samsvarer med uavhengighet

# Histogram over gjennomsnittene med tetthet for gjennomsnittene (N(0,1/sqrt(10))) 
hist(simres[,1],nclass=20,xlab="gj.snitt",prob=T,main="Histogram gj.snitt")
x=seq(-1.2,1.2,0.01)
lines(x,dnorm(x,0,1/sqrt(10)))
text(0.7,1.25,"X ~ N(0,1)")

# Histogram over empiriske varianser
hist(simres[,2],xlab="S^2",prob=T,main="Histogram emp.var",nclass=20)
x=seq(0,3,0.01)
lines(x,dgamma(x,9/2,9/2))


#Simulerer 1000 gammafordelte datasett med n=10 og ser på samsvar mellom gjennomsnitt og empirisk varians
simres=numeric(0)
for (sim in 1:1000){
  x=rgamma(10,3,2)
  simres=rbind(simres,c(mean(x),var(x)))
}
plot(simres[,1],simres[,2],xlab="gjennomsnitt",ylab="empirisk varians",main="Gammafordelte data")
cor(simres[,1],simres[,2])
# Med gammafordelte data har vi altså avhengighet mellom gjennomsnitt og empirisk varians

# Histogram over gjennomsnitt med gamma(3,0.5)-data
hist(simres[,1],xlab="gj.snitt",prob=T,main="Histogram over gj.snitt",nclass=20)
x=seq(0,3,0.01)
lines(x,dnorm(x,1.5,sqrt(0.75/10)))

# Histogram over empirisk varians med gamma(3,0.5)-data
hist(simres[,2],xlab="S^2",nclass=20,prob=T,main="Histogram emp.var")
x=seq(0,3,0.01)
lines(x,dgamma(x,9/2,9/2))
text(3,1,"Gammafordelte X")
# Gj.snittet er fortsatt normfordelt, men empirisk varians er ikke gamma (eller knyttet til kji-kvadrat)

## Tipsdata

# Leser inn data fra eksempel 8.17
tip = c(22.7,16.3,13.6,16.8,29.9,15.9,14.0,15.0,14.1,18.1,22.8,27.6,16.4,16.1,19.0,
      13.5,18.9,20.2,19.7,18.2,15.4,15.7,19.0,11.5,18.4,16.0,16.9,12.0,40.1,19.2)
#tip = tip[1:15]
n = length(tip)
mean.tip = mean(tip)
sd.tip = sd(tip)
c(mean.tip,sd.tip,n)

boxplot(tip)
# Histogram og normalt Q-Q-plott av dataene
#par(mfrow=c(1,2))
hist(tip)
abline(v=mean.tip,lty=2,lwd=2)
abline(v=median(tip),lty=3,lwd=2)
qqnorm(tip)
qqline(tip)


# Konfidensintervall basert paa normalfordeling
low.norm = mean.tip-qt(0.975,n-1)*sd.tip/sqrt(n)
up.norm = mean.tip+qt(0.975,n-1)*sd.tip/sqrt(n)
print(c(low.norm,up.norm),digits=3)

# Beregner bootstrap-estimat
B=1000
mean.boot=rep(0,B)
for (b in 1:B)
{
  index = sample(1:n,n,replace=TRUE)
  mean.boot[b] = mean(tip[index])
  #Alternativ kode
  #x.star = sample(tip,n,replace=TRUE)
  #mean.boot[b] = mean(x.star)
}

# Lager histogram og normalt Q-Q-plott av bootstrap-estimatene
# par(mfrow=c(1,2))
hist(mean.boot)
abline(v=mean(mean.boot),lty=2,lwd=2)
abline(v=median(mean.boot,lty=2,lwd=2))
qqnorm(mean.boot)
qqline(mean.boot)

# Beregner gjennomsnitt og standardavvik for bootstrap-estimatene
c(mean(mean.boot),sd(mean.boot))

# Beregner 95% konfidensinterval basert paa normalfordeling av bootstrap-estimatene
low.boot = mean(tip)-qnorm(0.975)*sd(mean.boot)
up.boot = mean(tip)+qnorm(0.975)*sd(mean.boot)
print(c(low.boot,up.boot),digits=3)

# Beregner 95% konfidensintervall etter persentilmetoden
mean.boot.sort=sort(mean.boot)
low.pers=mean.boot.sort[0.025*B]
up.pers=mean.boot.sort[0.975*B]
print(c(low.pers,up.pers),digits=3)

# R-pakken "boot" har prosedyrer for bootstrap-beregninger, deriblant BCa-intervallet
#install.packages("boot")
library(boot)

# Trekker bootstrap-utvalg. Merk at vi foerst maa definere en funksjon 
# som brukes til aa beregne bootstrap-estimatene (her gjennomsnitt)
my.mean = function(x, indices) 
{ 
  mean(x[indices])
}
mean.boot2= boot(tip, my.mean, R=1000)

# "mean.boot2" er et objekt med mange attibutter; disse kan en se ved aa skrive
names(mean.boot2)

# For eksempel er "mean.boot2$t0" estimatet basert på de 
# opprinnelige dataene, mens "mean.boot2$t" gir bootstrap-estimatene
mean.boot2$t0
summary(mean.boot2$t)

# BCa-konfidensintervallet:
boot.ci(mean.boot2,type=c("perc", "bca"))


# Bootstrap for medianer
my.median = function(x, indices) 
{ 
  median(x[indices])
}
median.boot2= boot(tip, my.median, R=1000)
boot.ci(median.boot2,type=c("perc", "bca"))


# KI binomisk p
# Eks Y = 5, n=20
y=30;n=100;phat=y/n
phat+c(-1,1)*1.96*sqrt(phat*(1-phat)/n)
prop.test(y,n)$conf.int
prop.test(y,n,correct=F)$conf.int

y=5;n=20;phat=y/n
phat+c(-1,1)*1.96*sqrt(phat*(1-phat)/n)
prop.test(y,n)$conf.int
prop.test(y,n,correct=F)$conf.int