data=c(rep(0,4),rep(1,6),rep(2,6),rep(3,8),rep(4,3),6,6,10)

#apriorni hustota - rovnomerna

verohodnost=function(x){
a=1
for (i in 1:length(data)){
a=a*dpois(data[i],x)
} 
return(a)
}

apost=function(x){
verohodnost(x)}

K=integrate(apost,0,Inf)$value

apost2=function(x){
apost2=apost(x)/K
return(apost2)}

curve(apost2,1,4,ylim=c(0,1.5),xlab=expression(lambda),
ylab=expression(paste(pi,"(",lambda,"|x)" )),main='Graf aposteriorni hustoty')

curve(dgamma(x,shape=77,scale=1/30),add=TRUE,col=7) #skutecna aposteriorni hustota
#1.)
#a)
f=function(x){
f=x*apost2(x)
}
integrate(f,0,Inf)$val[1]  #aposteriorni stredni hodnota

F50=function(x){
F50=integrate(apost2,0,x)$value[1]-0.5
return(F50)}

uniroot(F50,c(1,3))$root  #aposteriorni median

M=optimize(apost2,c(1,3),maximum=TRUE)$max  #aposteriorni modus
M

#################################
#b)

F25=function(x){
F25=integrate(apost2,0,x)$value[1]-0.025
return(F25)}

a=uniroot(F25,c(1,3))$root[1]
qgamma(0.025,shape=77,scale=1/30)

F975=function(x){
F975=integrate(apost2,0,x)$value[1]-0.975
return(F975)}

b=uniroot(F975,c(1,4))$root[1]
qgamma(0.975,shape=77,scale=1/30)

c(a,b)  # symetricky 95% verohodnostni interval

###########################################
#c)

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

F=function(x){
F=integrate(apost2,0,x)$value[1]
return(F)}


spolehlivost=function(c){

f=function(x){ #pomocna funkce pro hledani HPD intervalu pro dane c
f=apost2(x)-c
}
a=uniroot(f,c(1,M))$root
b=uniroot(f,c(M,4))$root

spolehlivost=F(b)-F(a)-0.95
return(spolehlivost)}

C=uniroot(spolehlivost, c(0.001,1))$root #vysledna hodnota c

f=function(x){
f=apost2(x)-C
}
a=uniroot(f,c(1,M))$root
b=uniroot(f,c(M,4))$root

c(a,b) # 95% highest posterior density interval



########################## predikce
#d)

predikce=function(y){
a=function(x){
a=apost2(x)*dpois(y,x)
}
val=integrate(a,0,Inf)[1]$value
return(val)
}

pred=rep(0,11)
ind=0:10
for (i in 0:10){
pred[i+1]=predikce(i)}

plot(pred~ind,type='h',xlab='x',ylab='p(x)',main='prediktivni pravd. funkce')


#2.)
opak=1000000
A=runif(opak,0,5)
B=runif(opak,0,2)

D=A[apost2(A)>B]  #realizace n.v. s aposteriornim rozdelenim
plot(density(D)) #jadrovy odhad hustoty nagenerovanych dat (apost. rozdeleni)
curve(apost2,col=2,add=TRUE)  #srovnani s teoretickou hustotou

mean(D)
median(D)

#jiny postup, jak odhadnout apost. stredni hodnotu je metoda MC z prednasky
X=runif(100000,0,1000)  # nahodny vyber z apriorniho rozdeleni
mean(X*verohodnost(X))/mean(verohodnost(X)) #odhad aposteriorni stredni hodnoty

#3.)

#a) \int_{0}^{1} x^9 dx = 1/10

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

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

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