# KOMMANDOER TIL FORELESNINGENE FOR UKE 10 

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

# MENINGSMAALING FRA JANUAR 2016
# n=721 har oppgitt hvilket parti de ville ha stemt paa
# x=214 ville ha stemt Arbeiderpariet

# Beregner likelihooden:
n=721
x=214
p=seq(0.25,0.35,0.001)
likelihood=choose(n,x)*p^x*(1-p)^(n-x)

# Plotter likelihood og log-likelihood
plot(p,likelihood,type='l')
plot(p,log(likelihood),type='l')

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


# EKSEMEL 1 FRA HEFTET TIL STORVIK: NEDBRYTING AV MYONER

# Leser inn dataene:

x = c(0.41040018,0.91061564,-0.61106896,0.39736684,0.37997637,0.34565436 ,
0.01906680,-0.28765977,-0.33169289,0.99989810,-0.35203164,0.10360470 ,
0.30573300,0.75283842,-0.33736278,-0.91455101,-0.76222116,0.27150040 ,
-0.01257456,0.68492778,-0.72343908,0.45530570,0.86249107,0.52578673 ,
0.14145264,0.76645754,-0.65536275,0.12497668,0.74971197,0.53839119)

# Histogram av dataene:
hist(x)

# Beregner likelihooden:
alpha=seq(-1,1,0.01)
likelihood=rep(0,length(alpha))
for (i in 1:length(alpha))
likelihood[i]= prod((1+alpha[i]*x)/2)

# Plotter likelihood og log-likelihood:
plot(alpha,likelihood,type='l')
plot(alpha,log(likelihood),type='l')




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

# EKSEMEL 12.14 I D&B: CHALLENGER ULYKKEN

# Leser inn dataene:
exmp12.14=read.table("http://www.uio.no/studier/emner/matnat/math/STK2120/v16/exmp12-14.txt",header=T)

# Trekker 70 fra alle temperaturm?lingene (som er i Farenheit-grader):
x=exmp12.14$temp-70
y=exmp12.14$fail

# Beregner og loglikelihood og likelihood:
b0 = seq(-3,0,0.1)
b1 = seq(-0.6,0,0.01)
loglik = matrix(nrow=length(b0),ncol=length(b1))
for(i in 1:length(b0))
for(j in 1:length(b1))
loglik[i,j] = sum((b0[i]+b1[j]*x)*y-log(1+exp(b0[i]+b1[j]*x)))
lik=exp(loglik)

#Plotter likelihooden
contour(b0,b1,lik)
persp(b0,b1,lik,,theta=330,phi=45,shade=1)

#Plotter likelihooden
contour(b0,b1,log(lik))
persp(b0,b1,log(lik),,theta=330,phi=45,shade=1)



# Plotter estimert sannsynligheten for feil som funksjon av temperatur
# (bruker at estimatene for beta0 og beta1 blir -1.2085 og -0.2322)
b0=-1.2085
b1=-0.2322
temp=seq(50,80,1)
plot(temp,exp(b0+b1*(temp-70))/(1+exp(b0+b1*(temp-70))),type='l', ylab="Sannsynlighet for feil")