n=20  #pocet hodu
n1=5  #pocet licu

verohodnost=function(theta){
a=theta^n1*(1-theta)^(n-n1)
return(a)}

curve(verohodnost,0,1,xlab=expression(theta),ylab="verohodnsot",  main="Verohodnost")
abline(v=n1/n,col=2)

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)

# bodove charakteristiky:

(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


# 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
krok=0.001
C=dbeta(M,alpha+n1,beta+n-n1)-krok
spolehlivost=0
while (spolehlivost<0.95){

f=function(x){
f=dbeta(x,alpha+n1, beta+n-n1)-C
}
x1=uniroot(f,c(0,M))$root
x2=uniroot(f,c(M,1))$root
spolehlivost=pbeta(x2,alpha+n1,beta+n-n1)-pbeta(x1,alpha+n-n1, beta+n-n1)
C=C-krok
}
spolehlivost
c(x1,x2) # highest posterior density interval

qbeta(0.975,alpha+n1,beta+n-n1)-qbeta(0.025,alpha+n1,beta+n-n1) #sirka equal-tail intervalu
x2-x1 #sirka HDP intervalu


# predikce budouciho pozorovani:

p0=(alpha+n-n1)/(alpha+beta+n)
p1=(alpha+n1)/(alpha+beta+n)
c(p0,p1)
plot(c(p0,p1)~c(0,1),ylim=c(0,1),type="h", xlab=" ", ylab="pravdepodobnost")
abline(h=0.5,col=3)


#####################
# vypocet integralu #
#####################
X=runif(1000000)
mean(X^9)


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


###############
# MC simulace #
###############
theta=rbeta(100000,alpha, beta)

In=mean(verohodnost(theta))
Jn=mean(verohodnost(theta)*theta)

Jn/In

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