#Monte Carlo estimation examples
x=seq(0,1,by=.01)
f=function(x) return(sin(x)*log(x+1)/x)
plot(x,f(x),type="l")

plot(x,f(x),type="l",ylim=c(0,1))

x=runif(100000)
y=runif(100000)
z=f(x)
mean(y<z)

#gamma mixture of Poissons (posterior predictive for Poisson with gamma prior)
lambda=rgamma(100000,1191+144/25,56+12/25) 
mean(lambda)
(1191+144/25)/(56+12/25)
var(lambda)
(1191+144/25)/(56+12/25)^2

y=rpois(100000,lambda)
mean(y)
var(y)
(1191+144/25)/(56+12/25)*(1+1/56)

median(y)
(sort(y))[50000]
(sort(y))[2500]
(sort(y))[97500]

#compare to Poisson
qpois(.025,(1191+144/25)/(56+12/25))
qpois(.975,(1191+144/25)/(56+12/25))

#Inverse CDF sampling pictoral demonstration
x=seq(-3.5,3.5,by=.01)
plot(x,pnorm(x),type="l")
u=runif(100)
for(i in 1:100) {
  points(0,u[i],pch=18)
  points(qnorm(u[i]),0,pch=17)
}
hist(u)
hist(qnorm(u))

u=runif(10000)
hist(u)
hist(qnorm(u))