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){
verohodnost<-prod(10*dpois(data,x))}

apost<-function(x){
verohodnost(x)*1
}

K<-integrate(Vectorize(apost),0,Inf)$value
apost2<-function(x){
apost2<-apost(x)/K
return(apost2)}

integrate(Vectorize(apost2),0,Inf) # kontrola, ze je to skutecne hustota

ind<-seq(1,4, by=0.01) #graf vykreslime na mrizce od 1 do 4
plot(Vectorize(apost2)(ind)~ind, type="l",xlab=expression(lambda),
ylab=expression(paste(pi,"(",lambda,"|x)" )),main='Graf aposteriorni hustoty',
ylim=c(0,1.5))

# je to hustota gamma rozdeleni s paramtry shape=77 a scale=1/30

#1.)
#a)
f<-function(x){
f<-x*apost2(x)
}
integrate(Vectorize(f),0,Inf)$val[1]  #aposteriorni stredni hodnota

F50<-function(x){
F50<-integrate(Vectorize(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(Vectorize(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(Vectorize(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(Vectorize(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){
return(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(Vectorize(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[Vectorize(apost2)(A)>B]  #realizace n.v. s aposteriornim rozdelenim
plot(density(D)) #jadrovy odhad hustoty nagenerovanych dat (apost. rozdeleni)
curve(dgamma(x,shape=77,scale=1/30),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 'temer' apriorniho rozdeleni
mean(X*Vectorize(verohodnost)(X))/mean(Vectorize(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))