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))
curve(dbeta(x,alpha+n1, beta+n-n1),0,1, col=2, xlab=expression(theta), main="Apriorni a aposteriorni rozdeleni", 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+n-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)


#####################
# vypocet integralu #
#####################
#\int_{0}^{1} x^9 dx = 1/10

X=runif(1000000)
mean(X^9)

#\int_{-\infty}^{\infty} e^{-|x|} dx = 2

X=rnorm(1000000)
f=function(x){
f=exp(-abs(x))/dnorm(x)
}
mean(f(X))


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