#--------------------------------------------------------------------
# Example of importance sampling:
#
#  We are interested in the mean value of exp(-abs(x)/tau)
#  with x ~ N(mu.p, sigma.p)
#
#  GdA 8/4/2016
#---------------------------------------------------------------------

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

n <- 100000 # number of events
s = 1000    # step between intermediate checks

# parameter of "f(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 <- 2

# 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()

# normal sampling
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
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')