#--------------------------------------------------------------------
# Example of importance sampling:
#
#  given X ~ N(mu.p, sigma.p) 
#  and a function of X, g(x) = exp(-abs(x)/tau),
#  we are interested in estimating E[g(x)]
#
#  (This example is only to show the difficulties   
#   using 'normal' Monte Carlo. The method becomes  
#   important in complicated cases!)
#  
#  GdA, May 2020
#--------------------------------------------------------------------

# function pause  -------------------------------------------
pause <- function() { cat ("\n >> Press enter to continue\n"); scan() }
#------------------------------------------------------------

n <- 10000 # number of events
s <- 100    # steps between intermediate checks

# parameter of "g(x)"  
tau=0.2

# parameters of the Gaussian "p(x)"
mu.p <- 3
sigma.p <- 1

# parameters of the Gaussian "q(x)"
mu.q <- 0
sigma.q <- 0.5

# plots
x <- seq(mu.p - 5*sigma.p, mu.p + 4*sigma.p, len=401)
plot(x, dnorm(x, mu.p, sigma.p), ty='l', ylim=range(exp(-abs(x)/tau)))
points(x, exp(-abs(x)/tau), ty='l', col='blue')
points(x,  dnorm(x, mu.q, sigma.q), ty='l', col='red')

pause()

# ordinary sampling, according to p(x)
x <- rnorm(n, mu.p, sigma.p)
fx <- exp(-abs(x)/tau)
ni <- seq(s, n, s)
Ifx = rep(0, length(ni))
for (i in 1:length(ni)) {
  Ifx[i] = sum(fx[1:ni[i]])/ni[i]
}

# importance sampling, according to q(x)
x <- rnorm(n, mu.q, sigma.q)
fx <- exp(-abs(x)/tau)*dnorm(x, mu.p, sigma.p)/dnorm(x, mu.q, sigma.q)
Ifx.IS = rep(0, length(ni))
for (i in 1:length(ni)) {
  Ifx.IS[i] = sum(fx[1:ni[i]])/ni[i]
}

plot(ni,Ifx, ylim=c(min(Ifx,Ifx.IS),max(Ifx,Ifx.IS)), col='black', ty='l')
points(ni,Ifx.IS, ty='l', col='red')