Ggplot2: гистограмма с нормальной кривой

Ответ 1

Думаю, я понял:

set.seed(1)
df <- data.frame(PF = 10*rnorm(1000))
ggplot(df, aes(x = PF)) + 
    geom_histogram(aes(y =..density..),
                   breaks = seq(-50, 50, by = 10), 
                   colour = "black", 
                   fill = "white") +
stat_function(fun = dnorm, args = list(mean = mean(df$PF), sd = sd(df$PF)))

Ответ 2

На этот вопрос ответили здесь и частично здесь.

Площадь под кривой плотности равна 1, а площадь под гистограммой равна ширине столбцов, умноженной на их высоту, т.е. ширина полосы умножается на общее количество не пропущенных наблюдений. Чтобы разместить оба элемента на одном и том же графике, необходимо изменить масштаб одного или другого, чтобы их области совпадали.

Если вы хотите, чтобы ось Y имела частоту отсчетов, есть несколько вариантов:

Сначала смоделируйте некоторые данные.

library(ggplot2)

set.seed(1)
dat_hist <- data.frame(
  group = c(rep("A", 200), rep("B",150)),
  value = c(rnorm(200, 20, 5), rnorm(150,25,10)))

# Set desired binwidth and number of non-missing obs
bw = 2
n_obs = sum(!is.na(dat_hist$value))

Вариант 1. Отобразите как гистограмму, так и кривую плотности как плотность, а затем измените масштаб оси y.

Это, пожалуй, самый простой подход для одной гистограммы. Используя подход, предложенный Карлосом, постройте гистограмму и кривую плотности как плотность

g <- ggplot(dat_hist, aes(value))  + 
geom_histogram(aes(y = ..density..), binwidth = bw, colour = "black") + 
stat_function(fun = dnorm, args = list(mean = mean(dat_hist$value), sd = sd(dat_hist$value)))

А затем измените масштаб оси Y.

ybreaks = seq(0,50,5) 
## On primary axis
g + scale_y_continuous("Counts", breaks = round(ybreaks / (bw * n_obs),3), labels = ybreaks)

## Or on secondary axis
g + scale_y_continuous("Density", sec.axis = sec_axis(
  trans = ~ . * bw * n_obs, name = "Counts", breaks = ybreaks))

Вариант 2. Масштабируйте кривую плотности с помощью функции stat_function

Код приведен в соответствие с ответом PatrickT.

ggplot(dat_hist, aes(value))  + 
  geom_histogram(colour = "black", binwidth = bw) + 
  stat_function(fun = function(x) 
    dnorm(x, mean = mean(dat_hist$value), sd = sd(dat_hist$value)) * bw * n_obs)

Вариант 3. Создание внешнего набора данных и графика с использованием geom_line.

В отличие от вышеуказанных опций, этот работает с фасетами. (ИЗМЕНЕНО для предоставления решения dplyr, а не plyr). Обратите внимание, что суммарный набор данных используется в качестве основного, а необработанные данные передаются только для гистограммы.

library(tidyverse)

dat_hist %>% 
  group_by(group) %>% 
  nest(value) %>% 
  mutate(y = map(data, ~ dnorm(
    .$value, mean = mean(.$value), sd = sd(.$value)
    ) * bw * sum(!is.na(.$value)))) %>% 
  unnest(data,y) %>% 

  ggplot(aes(x = value)) +
  geom_histogram(data = dat_hist, binwidth = bw, colour = "black") +
  geom_line(aes(y = y)) + 
  facet_wrap(~ group)

Вариант 4. Создание внешних функций для редактирования данных на лету

Может быть, немного сверх того, но может быть кому-то полезно?

## Function to create scaled dnorm data along full x axis range
dnorm_scaled <- function(data, x = NULL, binwidth = 1, xlim = NULL) {
  .x <- na.omit(data[,x])
  if(is.null(xlim))
    xlim = c(min(.x), max(.x))
  x_range = seq(xlim[1], xlim[2], length.out = 101)
  setNames(
    data.frame(
    x = x_range,
    y = dnorm(x_range, mean = mean(.x), sd = sd(.x)) * length(.x) * binwidth),
    c(x, "y"))
}

## Function to apply over groups
dnorm_scaled_group <- function(data, x = NULL, group = NULL, binwidth = NULL, xlim = NULL) {
  dat_hists <- lapply(
    split(data, data[, group]), dnorm_scaled,
      x = x, binwidth = binwidth, xlim = xlim)
  for(g in names(dat_hists))
    dat_hists[[g]][, "group"] <- g
  setNames(do.call(rbind, dat_hists), c(x, "y", group))
}

## Single histogram
ggplot(dat_hist, aes(value)) + 
  geom_histogram(binwidth = bw, colour = "black") + 
  geom_line(data = ~ dnorm_scaled(., "value", binwidth = bw), 
            aes(y = y)) 

## With a single faceting variable
ggplot(dat_hist, aes(value))  + 
  geom_histogram(binwidth = 2, colour = "black") + 
  geom_line(data = ~ dnorm_scaled_group(
    ., x = "value", group = "group", binwidth = 2, xlim = c(0,50)), 
    aes(y = y)) +
  facet_wrap(~ group)

Ответ 3

Это расширенный комментарий к ответу Дж. Виллимана. Я нашел ответ J очень полезным. Во время игры я обнаружил способ упростить код. Я не говорю, что это лучший способ, но я думал, что упомяну это.

Обратите внимание, что ответ JWilliman предоставляет счет на оси Y и "хак" для масштабирования соответствующего приближения нормальной плотности (которое в противном случае охватило бы общую площадь 1 и, следовательно, имело бы намного более низкий пик).

Основной смысл этого комментария: упрощенный синтаксис внутри stat_function, путем передачи необходимых параметров в функцию эстетики, например,

aes(x = x, mean = 0, sd = 1, binwidth = 0.3, n = 1000)

Это позволяет избежать необходимости передавать args = в stat_function и, следовательно, более удобно для пользователя. Ладно, это не сильно отличается, но, надеюсь, кто-то найдет это интересным.

# parameters that will be passed to ''stat_function''
n = 1000
mean = 0
sd = 1
binwidth = 0.3 # passed to geom_histogram and stat_function
set.seed(1)
df <- data.frame(x = rnorm(n, mean, sd))

ggplot(df, aes(x = x, mean = mean, sd = sd, binwidth = binwidth, n = n)) +
    theme_bw() +
    geom_histogram(binwidth = binwidth, 
        colour = "white", fill = "cornflowerblue", size = 0.1) +
stat_function(fun = function(x) dnorm(x, mean = mean, sd = sd) * n * binwidth,
    color = "darkred", size = 1)

Ответ 4

Этот код должен сделать это:

set.seed(1)
z <- rnorm(1000)

qplot(z, geom = "blank") + 
geom_histogram(aes(y = ..density..)) + 
stat_density(geom = "line", aes(colour = "bla")) + 
stat_function(fun = dnorm, aes(x = z, colour = "blabla")) + 
scale_colour_manual(name = "", values = c("red", "green"), 
                               breaks = c("bla", "blabla"), 
                               labels = c("kernel_est", "norm_curv")) + 
theme(legend.position = "bottom", legend.direction = "horizontal")

Примечание. Я использовал qplot, но вы можете использовать более универсальный ggplot.