27

I wrote the following code in R 

x=seq(-7,10,length=200)
y1=dnorm(x,mean=0,sd=1)
plot(x,y1,type="l",lwd=2,col="red")
y2=dnorm(x,mean=3,sd=2)
lines(x,y2,type="l",lwd=2,col="blue")

How can I shade the area under both curves (known as the overlap between the two curves).

I will highly appreciate any suggestions.

Andrie
  • 176,377
  • 47
  • 447
  • 496
Muna
  • 271
  • 1
  • 3
  • 3

3 Answers3

30

Oh, well, @SachaEpskamp beat me to it, but here is my much less elegant solution.

shade_under_curve <- function(fun, xmin, xmax, length=100){
  xvals <- seq(xmin, xmax, length=length)
  dvals <- match.fun(fun)(xvals)
  polygon(c(xvals,rev(xvals)),c(rep(0,length),rev(dvals)),col="gray")
}


y1 <- function(x)sapply(x, function(xt)dnorm(xt,mean=0,sd=1))
y2 <- function(x)sapply(x, function(xt)dnorm(xt,mean=3,sd=2))

my.fun <- function(x){sapply(x, function(xt)min(y1(xt), y2(xt)))}

Edit to include initial plot:

plot(y1, -10, 10, col="red")
curve(y2, add=TRUE, col="blue")
shade_under_curve(my.fun, -10, 10, length=1000)

enter image description here

Andrie
  • 176,377
  • 47
  • 447
  • 496
  • Andrie thank you so much for the help, but when I tried to run the code I couldn't get the whole graph as you have it up here – Muna Jul 22 '11 at 08:53
  • @Amalhelu, Thank you for pointing this out. I have edited the plot code. It should now work. (I had forgotten that the `curve` doesn't create a new plot). – Andrie Jul 22 '11 at 09:04
25

Add the following line:

polygon(x,pmin(y1,y2),col="gray")

This basically works exactly like the pen tool in photoshop, where the first vector, x are the x-coordinates and the second vector, pmin(y1,y2) are the y coordinates. pmin gives you a vector with the minimal values of two vectors elementwise, which corresponds to the y coordinates of the top of the overlap.

EDIT:

I prefer using curve() (as Andrie suggested), which can be used to plot a function. You can save it's coordinates while plotting too and us it in exactly the same way:

fun1 <- curve(dnorm(x,mean=0,sd=1),type="l",lwd=2,col="red")
fun2 <- curve(dnorm(x,mean=3,sd=2),type="l",lwd=2,col="blue",add=TRUE)
polygon(fun1$x,pmin(fun1$y,fun2$y),col="gray")
Sacha Epskamp
  • 46,463
  • 20
  • 113
  • 131
11

Here is a solution using ggplot2

library(ggplot2)
x  = seq(-7, 10, length = 200)
y1 = dnorm(x, mean = 0,sd = 1)
y2 = dnorm(x, mean = 3,sd = 2)

mydf = data.frame(x, y1, y2)

p0 = ggplot(mydf, aes(x = x)) +                         
     geom_line(aes(y = y1), colour = 'blue') +
     geom_line(aes(y = y2), colour = 'red') +
     geom_area(aes(y = pmin(y1, y2)), fill = 'gray60')
Ramnath
  • 54,439
  • 16
  • 125
  • 152
  • This is great, but if I reduce the number of x values to, for example, 20 and re-plot, the polygon no longer neatly fits the intersection of the two curves. I guess this is because its vertices are defined by the x-y values in the data, however the intersection of the two lines occurs between two x-y values. Any suggestions on what I should do in this situation? – Nicholas Oct 22 '15 at 11:42