## 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.3731862

Create a plot of the target and proposal distributions

df1 <- data.frame(x, q, g) %>%
pivot_longer(cols = !x, names_to = "grp", values_to = "p")
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('g(theta|y)', 'q(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)

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDEwLjIiCmF1dGhvcjogIkFraSBWZWh0YXJpLCBNYXJrdXMgUGFhc2luaWVtaSIKZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdGhlbWU6IHJlYWRhYmxlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCi0tLQojIyBJbXBvcnRhbmNlIHNhbXBsaW5nIGV4YW1wbGUKCmdncGxvdDIgYW5kIGdyaWRFeHRyYSBhcmUgdXNlZCBmb3IgcGxvdHRpbmcsIHRpZHlyIGZvciBtYW5pcHVsYXRpbmcKZGF0YSBmcmFtZXMKCmBgYHtyIHNldHVwLCBtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShnZ3Bsb3QyKQp0aGVtZV9zZXQodGhlbWVfbWluaW1hbCgpKQpsaWJyYXJ5KGdyaWRFeHRyYSkKbGlicmFyeSh0aWR5cikKYGBgCgpGYWtlIGludGVyZXN0aW5nIGRpc3RyaWJ1dGlvbgoKYGBge3IgfQp4IDwtIHNlcSgtNCwgNCwgbGVuZ3RoLm91dCA9IDIwMCkKciA8LSBjKDEuMSwgMS4zLCAtMC4xLCAtMC43LCAwLjIsIC0wLjQsIDAuMDYsIC0xLjcsCiAgICAgICAxLjcsIDAuMywgMC43LCAxLjYsIC0yLjA2LCAtMC43NCwgMC4yLCAwLjUpCmBgYAoKQ29tcHV0ZSB1bm5vcm1hbGl6ZWQgdGFyZ2V0IGRlbnNpdHkgKG5hbWVkIHEsIHRvIGVtcGhhc2l6ZSB0aGF0IGl0CmRvZXMgbm90IG5lZWQgdG8gYmUgbm9ybWFsaXplZCkuCgpgYGB7ciB9CnEgPC0gZGVuc2l0eShyLCBidyA9IDAuNSwgbiA9IDIwMCwgZnJvbSA9IC00LCB0byA9IDQpJHkKYGBgCgpHYXVzc2lhbiBwcm9wb3NhbCBkaXN0cmlidXRpb24KCmBgYHtyIH0KZyA8LSBkbm9ybSh4KQp3IDwtIHEvZwpycyA8LSBybm9ybSgxMDApCiMgZmluZCBuZWFyZXN0IHBvaW50IGZvciB3aGljaCB0aGUga2VybmVsIGhhcyBiZWVuIGV2YWx1YXRlZCBmb3IgZWFjaCBzYW1wbGUKcnNpIDwtIHNhcHBseShycywgZnVuY3Rpb24oYXJnKSB3aGljaC5taW4oYWJzKGFyZyAtIHgpKSkKYGBgCgpTZWxmLW5vcm1hbGl6ZWQgaW1wb3J0YW5jZSB3ZWlnaHRzIGFuZCB0aGUgZXhwZWN0YXRpb24gd3J0IHEKCmBgYHtyIH0Kd3IgPC0gcVtyc2ldL2Rub3JtKHhbcnNpXSkKd3JuIDwtIHdyL3N1bSh3cikKKEV4IDwtIHN1bSh3cm4qeFtyc2ldKSkKYGBgCgpDcmVhdGUgYSBwbG90IG9mIHRoZSB0YXJnZXQgYW5kIHByb3Bvc2FsIGRpc3RyaWJ1dGlvbnMKCmBgYHtyIH0KZGYxIDwtIGRhdGEuZnJhbWUoeCwgcSwgZykgJT4lCiAgcGl2b3RfbG9uZ2VyKGNvbHMgPSAheCwgbmFtZXNfdG8gPSAiZ3JwIiwgdmFsdWVzX3RvID0gInAiKQpkaXN0ciA8LSBnZ3Bsb3QoZGF0YSA9IGRmMSkgKwogIGdlb21fbGluZShhZXMoeCwgcCwgZmlsbCA9IGdycCwgY29sb3IgPSBncnApKSArCiAgbGFicyh0aXRsZSA9ICdUYXJnZXQgYW5kIHByb3Bvc2FsIGRpc3RyaWJ1dGlvbnMnLCB4ID0gJycsIHkgPSAnJykgKwogIHNjYWxlX2NvbG9yX2Rpc2NyZXRlKGxhYmVscyA9IGMoJ2codGhldGF8eSknLCAncSh0aGV0YSknKSkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICdib3R0b20nLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKQ3JlYXRlIGEgcGxvdCBvZiB0aGUgc2FtcGxlcyBhbmQgaW1wb3J0YW5jZSB3ZWlnaHRzCgpgYGB7ciB9CnNhbXAgPC0gZ2dwbG90KCkgKwogIGdlb21fbGluZShhZXMoeCwgdywgY29sb3IgPSAnMScpKSArCiAgZ2VvbV9zZWdtZW50KGFlcyh4ID0geFtyc2ldLCB4ZW5kID0geFtyc2ldLCB5ID0gMCwgeWVuZCA9IHdyKSwKICAgICAgICAgICAgICAgYWxwaGEgPSAwLjUsIGNvbG9yID0gJ3N0ZWVsYmx1ZScpICsKICBsYWJzKHRpdGxlID0gJ1NhbXBsZXMgYW5kIGltcG9ydGFuY2Ugd2VpZ2h0cycsIHggPSAnJywgeSA9ICcnKSArCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGMoJ3N0ZWVsYmx1ZScpLCBsYWJlbHMgPSAncSh0aGV0YXx5KS9nKHRoZXRhKScpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKCkNvbWJpbmUgdGhlIHBsb3RzCgpgYGB7ciB9CmdyaWQuYXJyYW5nZShkaXN0ciwgc2FtcCkKYGBgCgo=