## Importance sampling example

ggplot2 and gridExtra are used for plotting, tidyr for manipulating data frames

``````library(ggplot2)
theme_set(theme_minimal())
library(gridExtra)
library(tidyr)``````

Fake interesting distribution

``````x <- seq(-4, 4, length.out = 200)
r <- c(1.1, 1.3, -0.1, -0.7, 0.2, -0.4, 0.06, -1.7,
1.7, 0.3, 0.7, 1.6, -2.06, -0.74, 0.2, 0.5)``````

Compute unnormalized target density (named q, to emphasize that it does not need to be normalized).

``q <- density(r, bw = 0.5, n = 200, from = -4, to = 4)\$y``

Gaussian proposal distribution

``````g <- dnorm(x)
w <- q/g
rs <- rnorm(100)
# find nearest point for which the kernel has been evaluated for each sample
rsi <- sapply(rs, function(arg) which.min(abs(arg - x)))``````

Self-normalized importance weights and the expectation wrt q

``````wr <- q[rsi]/dnorm(x[rsi])
wrn <- wr/sum(wr)
(Ex <- sum(wrn*x[rsi]))``````
``## [1] 0.09706821``

Create a plot of the target and proposal distributions

``````df1 <- data.frame(x, q, g) %>% gather(grp, p, -x)
distr <- ggplot(data = df1) +
geom_line(aes(x, p, fill = grp, color = grp)) +
labs(title = 'Target and proposal distributions', x = '', y = '') +
scale_color_discrete(labels = c('q(theta|y)', 'g(theta)')) +
theme(legend.position = 'bottom', legend.title = element_blank())``````
``## Warning: Ignoring unknown aesthetics: fill``

Create a plot of the samples and importance weights

``````samp <- ggplot() +
geom_line(aes(x, w, color = '1')) +
geom_segment(aes(x = x[rsi], xend = x[rsi], y = 0, yend = wr),
alpha = 0.5, color = 'steelblue') +
labs(title = 'Samples and importance weights', x = '', y = '') +
scale_color_manual(values = c('steelblue'), labels = 'q(theta|y)/g(theta)') +
theme(legend.position = 'bottom', legend.title = element_blank())``````

Combine the plots

``grid.arrange(distr, samp)``

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDEwLjIiCmF1dGhvcjogIkFraSBWZWh0YXJpLCBNYXJrdXMgUGFhc2luaWVtaSIKZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdGhlbWU6IHJlYWRhYmxlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCi0tLQojIyBJbXBvcnRhbmNlIHNhbXBsaW5nIGV4YW1wbGUKCmdncGxvdDIgYW5kIGdyaWRFeHRyYSBhcmUgdXNlZCBmb3IgcGxvdHRpbmcsIHRpZHlyIGZvciBtYW5pcHVsYXRpbmcKZGF0YSBmcmFtZXMKCmBgYHtyIHNldHVwLCBtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFy