n<-20  #pocet hodu
n1<-5  #pocet licu

#definujeme si verohodnostni funkci
verohodnost<-function(theta){
a<-theta^n1*(1-theta)^(n-n1)
return(a)}

curve(verohodnost,0,1,xlab=expression(theta),ylab="verohodnsot",  main="Verohodnost") #graf verohodnosti
abline(v=n1/n,col=2) #pridame bod maxima (maximalne verohodny odhad)

alpha<-100  #parametry apriorniho beta rozdeleni
beta<-100   #parametry apriorniho beta rozdeleni

curve(dbeta(x,alpha, beta),0,1, ylim=c(0,12), xlab=expression(theta), main="Apriorni a aposteriorni rozdeleni")
curve(dbeta(x,alpha+n1, beta+n-n1),0,1, col=2, xlab=expression(theta), add=T)
legend("topright", legend=c("apriorni","aposteriorni"),col=1:2,lty=1)

# bodove odhady:

(alpha+n1)/(alpha+beta+n) # aposteriorni stredni hodnota
(alpha+n1-1/3)/(alpha+beta+n-2/3) # priblizny aposteriorni median
(alpha+n1-1)/(alpha+beta+n-2) # aposteriorni modus


# 95% verohodnostni intervaly:

c(qbeta(0.025,alpha+n1,beta+n-n1),qbeta(0.975,alpha+n1,beta+n-n1)) # equal-tail interval

M<-(alpha+n1-1)/(alpha+beta+n-2) #aposteriorni modus
H<-dbeta(M,alpha+n1,beta+n-n1)   #hodnota hustoty v aposteriornim modu

#definujeme si funkci, ktera pro danou hodnotu c spocita prislusnou spolehlivost-0.95

spolehlivost<-function(c){

f<-function(x){ #pomocna funkce pro hledani HPD intervalu pro dane c
f<-dbeta(x,alpha+n1, beta+n-n1)-c
}
a<-uniroot(f,c(0,M))$root
b<-uniroot(f,c(M,1))$root

spolehlivost<-pbeta(b,alpha+n1,beta+n-n1)-pbeta(a,alpha+n1, beta+n-n1)-0.95
return(spolehlivost)}

C=uniroot(spolehlivost, c(0,H))$root #vysledna hodnota c

f<-function(x){
f<-dbeta(x,alpha+n1, beta+n-n1)-C
}
a<-uniroot(f,c(0,M))$root
b<-uniroot(f,c(M,1))$root

c(a,b) # highest posterior density interval

qbeta(0.975,alpha+n1,beta+n-n1)-qbeta(0.025,alpha+n1,beta+n-n1) #sirka equal-tail intervalu
b-a #sirka HPD intervalu


# predikce budouciho pozorovani:

p0<-(alpha+n-n1)/(alpha+beta+n) #predikce rubu
p1<-(alpha+n1)/(alpha+beta+n)   #predikce lice
c(p0,p1)
plot(c(p0,p1)~c(0,1),ylim=c(0,1),type="h", xlab="x_{n+1}", ylab="pravdepodobnost")
abline(h=0.5,col=3)


###############
# MC simulace #
###############
#priklad na vypocet aposteriorni stedni hodnoty beta rozdeleni

theta<-rbeta(100000,alpha, beta)

In<-mean(verohodnost(theta))
Jn<-mean(verohodnost(theta)*theta)

Jn/In #aposteriorni stredni hodnota ziskana pomoci simulaci

f<-function(x){
f<-verohodnost(x)*dbeta(x,alpha, beta)
}
g<-function(x){
g<-x*verohodnost(x)*dbeta(x,alpha, beta)
}


integrate(g,0,1)$val/integrate(f,0,1)$val #aposteriorni stredni hodnota ziskana pomoci numerickeho integrovani