Plot two densities with shaded area under the curve

plot_two_densities.R

# http://www.neilmclatchie.com/using-r-shading-area-under-a-distribution/
stp <- 0.01  # "resolution"
left.limit <- -6
right.limit <- 6

plot(
  x = NA,
  type = "n",
  lwd = 2,
  lty = 1,
  col = "black",
  xlim = c(-6, 10),
  ylim = c(0, 0.4),
  frame.plot = FALSE,
  xlab = " ",
  ylab = " ",
  axes = FALSE,
  main = paste(" ")
)
Axis(side = 1, at = seq(-6, 10, by = 1))

# LEFT curve
c1.mean <- 0
c1.sd <- 1
start1 <- qnorm(0.95, mean = c1.mean, sd = c1.sd)
lx <- seq(from = start1, to = right.limit, by = stp)
ly <- dnorm(x = lx, mean = c1.mean, sd = c1.sd)
ly[1] <- 0  # very important to set this to zero to pin polygon to x axis

curve(dnorm(x, mean = c1.mean, sd = c1.sd), add = TRUE)
polygon(x = lx, y = ly, col = "grey80")

head(data.frame(lx, ly))
# RIGHT curve
c2.mean <- 3
c2.sd <- 1
stop2 <- qnorm(0.05, mean = c2.mean, sd = c2.sd)
rx <- seq(from = left.limit, to = stop2, by = stp)
ry <- dnorm(rx, mean = c2.mean, sd = c2.sd)
ry[length(ry)] <- 0  # ibidem

curve(dnorm(x, mean = c2.mean, sd = c2.sd), add = TRUE)
polygon(x = rx, y = ry, density = 10)

# Add vertical line
vline <- 1.5
vline.height <- 0.4
segments(x0 = vline, y0 = -0.1, x1 = vline, y1 = vline.height,
         col = "red")