Overlay histogram with density curve

Question

I am trying to make a histogram of density values and overlay that with the curve of a density function (not the density estimate).

Using a simple standard normal example, here is some data:

x <- rnorm(1000)

I can do:

q <- qplot( x, geom="histogram")
q + stat_function( fun = dnorm )

but this gives the scale of the histogram in frequencies and not densities. with ..density.. I can get the proper scale on the histogram:

q <- qplot( x,..density.., geom="histogram")
q

But now this gives an error:

q + stat_function( fun = dnorm )

Is there something I am not seeing?

Another question, is there a way to plot the curve of a function, like curve(), but then not as layer?

Driftage

The issue is that you have defined a global y for your plot using ..density.. inside qplot. This confuses stat_function. The easiest fix would be to write qplot(x, geom = 'blank') + geom_histogram(aes(y = ..density..)) + stat_function(fun = dnorm). See my detailed answer below – Flesh

The equivalent to curve(dnorm, -4, 4) would be qplot(x = -4:4, stat = 'function', fun = dnorm, geom = 'line') – Flesh

Ah right, I tried that with the function as first argument but see now what went wrong. Thanks! – Driftage

Answer 1

Here you go!

# create some data to work with
x = rnorm(1000);

# overlay histogram, empirical density and normal density
p0 = qplot(x, geom = 'blank') +   
  geom_line(aes(y = ..density.., colour = 'Empirical'), stat = 'density') +  
  stat_function(fun = dnorm, aes(colour = 'Normal')) +                       
  geom_histogram(aes(y = ..density..), alpha = 0.4) +                        
  scale_colour_manual(name = 'Density', values = c('red', 'blue')) + 
  theme(legend.position = c(0.85, 0.85))

print(p0)

Flesh

P.S. If one works with real data, make sure to pass the empirical mean and sd arguments to dnorm function, see stat_function help for syntax. – Skirting

Just out of curiosity: How would this be done using the ggplot() function? I just barely understood the way ggplot() works, so I feel a little weird using this approach for my stuff. – Solita

@Solita you could swap the first line out for something like this "ggplot(data.frame(x), aes(x=x)) +" – Donofrio

@Solita Why is that? Without passing mean and sd in args to stat_function I get nothing at all. – Sclar

There is a problem with overlaying histograms and density estimations, which is that the density estimations should really be shifted half a binwidth to make for the most accurate and aesthetically pleasing presentation. I have not been able to figure out how to do this. Any takers? – Morganite

Could you give line by line explanation please? How does ..density.. work? Copy your code to R works, but changing to my own data never works – Llywellyn

Answer 2

A more bare-bones alternative to Ramnath's answer, passing the observed mean and standard deviation, and using ggplot instead of qplot:

df <- data.frame(x = rnorm(1000, 2, 2))

# overlay histogram and normal density
ggplot(df, aes(x)) +
  geom_histogram(aes(y = after_stat(density))) +
  stat_function(
    fun = dnorm, 
    args = list(mean = mean(df$x), sd = sd(df$x)), 
    lwd = 2, 
    col = 'red'
  )

Lustful

This is a very convenient answer, as it provides a way to plot a histogram and a density curve even when they belong to different distributions, if needed (as it was for me). Thank you! – Multiply

The original question is about fitting a density curve, not specifically a single Gaussian. If you want to see why this solution doesn't work, try setting the data to df <- data.frame(x = c(rnorm(1000, 2, 2), rnorm(1000, 12, 2))) – Reed

@Megatron, No, OP asked for "the curve of a density function (not the density estimate)". So I still think this is correct. Your example shows that the normal density function may not be a good description of the data in some cases, but that is besides the point. – Lustful

Answer 3

What about using geom_density() from ggplot2? Like so:

df <- data.frame(x = rnorm(1000, 2, 2))

ggplot(df, aes(x)) +
  geom_histogram(aes(y=..density..)) +  # scale histogram y
  geom_density(col = "red")

This also works for multimodal distributions, for example:

df <- data.frame(x = c(rnorm(1000, 2, 2), rnorm(1000, 12, 2), rnorm(500, -8, 2)))

ggplot(df, aes(x)) +
  geom_histogram(aes(y=..density..)) +  # scale histogram y
  geom_density(col = "red")

Floriated

Because OP asked for "the curve of a density function (not the density estimate)". geom_density gives the density estimate. – Lustful

Maybe not what the OP asked for, but this did help with what I was looking for! – Seedbed

@Lustful What is the difference between density function and density estimate? – Dovecote

Answer 4

I'm trying for iris data set. You should be able to see graph you need in these simple code:

ker_graph <- ggplot(iris, aes(x = Sepal.Length)) + 
geom_histogram(aes(y = ..density..),
colour = 1, fill = "white") +
geom_density(lwd = 1.2,
linetype = 2,
colour = 2)

Misdeem