В ggplot, как заполнить область между двумя нормальными кривыми

У меня есть две нормальные кривые, и я хочу заполнить правую область между обеими кривыми, поэтому левая кривая является нижней границей y, а правая кривая — верхней границей y. Для построения кривых я использую stat_function(), поэтому ggplot рисует кривую без определения столбца y в aes(). Я нарисовал область заливки между кривой и осью X, но мне нужна область между обеими кривыми, и трюк с опустошением левой кривой с помощью NA, похоже, не работает, как я ожидал.

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

Как мне это сделать?

library(ggplot2)
library(ggthemes)

graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2) {

  Xmin1 <- Xmedia1-4*Xdt1
  Xmax1 <- Xmedia1+4*Xdt1

  Xmin2 <- Xmedia2-4*Xdt2
  Xmax2 <- Xmedia2+4*Xdt2

  Ymax1 <- max(dnorm(Xmedia1, Xmedia1, Xdt1))
  Ymax2 <- max(dnorm(Xmedia2, Xmedia2, Xdt2))

  Xmin <- min(Xmin1, Xmin2)
  Xmax <- max(Xmax1, Xmax2)

  ggplot(data.frame(X = c(Xmin, Xmax)), aes(x = X)) +
  geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
  stat_function(fun = dnorm, 
              args = c(Xmedia1, Xdt1), 
              linewidth = 1, 
              colour = "grey") +
  stat_function(fun = dnorm, 
              args = c(Xmedia2, Xdt2), 
              linewidth = 1, 
              colour = "black") +
  geom_segment(aes(x = Xmedia1, y = 0, xend = Xmedia1, yend = Ymax1), 
             linetype = "dashed", 
             linewidth = 0, 
             colour = "grey") +
  geom_segment(aes(x = Xmedia2, y = 0, xend = Xmedia2, yend = Ymax2), 
             linetype = "dashed", 
             linewidth = 0, 
             colour = "black") +
  ####################################################################
  stat_function(fun = dnorm,
            args = c(Xmedia2, Xdt2),
            xlim = c(Xmedia2+1.5*Xdt2,Xmax2),
            geom = "area",
            fill = "red",
            alpha = 0.5) +
  stat_function(fun = dnorm,
              args = c(Xmedia1, Xdt1),
              xlim = c(Xmedia1,Xmax1),
              geom = "area",
              fill = NA,
              alpha = 0.01) +

  ##################################################################

  theme(
    line = element_blank(),
    axis.line.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis. Ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend. Position = "none",
    panel. Grid = element_blank(),
    panel. Background = element_rect(fill = "lightgray", colour = NA),
 ) +
 xlim(c(Xmin, Xmax)) 
}

g1 <- graf_normal(250, 7, 253, 7)

g1

График обеих кривых, которые я получаю, таков:

Спасибо,

Обновлено: Используя код @stephan и играя с фильтрацией данных, я смог сделать это проще, используя geom_ribbon():

Классный способ различать перекрывающиеся зоны! Полный код:

library(ggplot2)

graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2, n = 1000) {
  x1 <- Xmedia1 + 4 * Xdt1 * seq(-1, 1, length. Out = n)
  x2 <- Xmedia2 + 4 * Xdt2 * seq(-1, 1, length. Out = n)

  dat <- data. Frame(
    x = union(x1, x2)
  )
  dat$y1 <- dnorm(dat$x, Xmedia1, Xdt1)
  dat$y2 <- dnorm(dat$x, Xmedia2, Xdt2)

  Ymax1 <- dnorm(Xmedia1, Xmedia1, Xdt1)
  Ymax2 <- dnorm(Xmedia2, Xmedia2, Xdt2)

  ggplot(dat, aes(x)) +
    geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
    geom_ribbon(
      data =  subset(dat, x >= Xmedia2 + 1.5 * Xdt2),
      aes(ymin = y1, ymax = y2),
      fill = "red", alpha = 0.8
    ) +
    geom_ribbon(
      data =  subset(dat, (x <= Xmedia2 + 1.5 * Xdt2) & (y2 > y1)),
      aes(ymin = y1, ymax = y2),
      fill = "red", alpha = 0.2
    ) +
    geom_ribbon(
      data = subset(dat, x <= Xmedia1 - 1.5 * Xdt2),
      aes(ymin = y1, ymax = y2),
      fill = "blue", alpha = 0.8
    ) +
    geom_ribbon(
      data = subset(dat, (x <= Xmedia2 ) & (y1 > y2)),
      aes(ymin = y1, ymax = y2),
      fill = "blue", alpha = 0.2
    ) +
    annotate(
      geom = "segment",
      x = c(Xmedia1, Xmedia2), y = 0,
      xend = c(Xmedia1, Xmedia2), yend = c(Ymax1, Ymax2),
      linetype = "dashed",
      linewidth = 1,
      colour = c("grey", "black")
    ) +
    geom_line(aes(y = y1), linewidth = 1, colour = "grey") +
    geom_line(aes(y = y2), linewidth = 1, colour = "black") +
    theme(
      line = element_blank(),
      axis.line.y = element_blank(),
      axis. Text = element_blank(),
      axis. Ticks = element_blank(),
      axis. Title = element_blank(),
      legend. Position = "none",
      panel. Grid = element_blank(),
      panel. Background = element_rect(fill = "lightgray", colour = NA),
    )
}

graf_normal(250, 7, 253, 7)

Однако код работает не для всех кривых, работающих над ним!:

graf_normal(250, 7, 253, 3)

r ggplot2

19.02.2023 19:02

Стоит ли изучать PHP в 2023-2024 годах?

Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...

Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией

В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.

Приемы CSS-макетирования - floats и Flexbox

Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...

Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest

В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...

Концепция локализации и ее применение в приложениях React ⚡️

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

Пользовательский скаляр GraphQL

Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...

Перейти к ответу Данный вопрос помечен как решенный

Ответы 2

Ответ принят как подходящий

Одним из вариантов заполнения области между кривыми нормалей может быть использование ggh4x::stat_difference, что, однако, требует вычисления значений плотности вручную и рисования с помощью geom_line вместо того, чтобы полагаться на stat_function():

library(ggplot2)
library(ggh4x)
graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2, n = 101) {
  x1 <- Xmedia1 + 4 * Xdt1 * seq(-1, 1, length.out = n)
  x2 <- Xmedia2 + 4 * Xdt2 * seq(-1, 1, length.out = n)
  
  dat <- data.frame(
    x = union(x1, x2)
  )
  dat$y1 <- dnorm(dat$x, Xmedia1, Xdt1)
  dat$y2 <- dnorm(dat$x, Xmedia2, Xdt2)

  Ymax1 <- dnorm(Xmedia1, Xmedia1, Xdt1)
  Ymax2 <- dnorm(Xmedia2, Xmedia2, Xdt2)

  ggplot(dat, aes(x)) +
    geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
    ggh4x::stat_difference(
      data = ~ subset(.x, x >= Xmedia2 + 1.5 * Xdt2),
      aes(ymin = y1, ymax = y2)
    ) +
    annotate(
      geom = "segment",
      x = c(Xmedia1, Xmedia2), y = 0,
      xend = c(Xmedia1, Xmedia2), yend = c(Ymax1, Ymax2),
      linetype = "dashed",
      linewidth = 1,
      colour = c("grey", "black")
    ) +
    geom_line(aes(y = y1), linewidth = 1, colour = "grey") +
    geom_line(aes(y = y2), linewidth = 1, colour = "black") +
    scale_fill_manual(values = c(scales::alpha("red", .5), "transparent")) +
    theme(
      line = element_blank(),
      axis.line.y = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      axis.title = element_blank(),
      legend.position = "none",
      panel.grid = element_blank(),
      panel.background = element_rect(fill = "lightgray", colour = NA),
    )
}

graf_normal(250, 7, 253, 7)

РЕДАКТИРОВАТЬ На самом деле stat_differnce не требуется для вашего случая. Думал слишком сложно. Как отметил @JuanRiera в своем комментарии, мы могли бы заполнить область с помощью geom_ribbon:

library(ggplot2)

graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2, n = 101) {
  x1 <- Xmedia1 + 4 * Xdt1 * seq(-1, 1, length.out = n)
  x2 <- Xmedia2 + 4 * Xdt2 * seq(-1, 1, length.out = n)

  dat <- data.frame(
    x = union(x1, x2)
  )
  dat$y1 <- dnorm(dat$x, Xmedia1, Xdt1)
  dat$y2 <- dnorm(dat$x, Xmedia2, Xdt2)

  Ymax1 <- dnorm(Xmedia1, Xmedia1, Xdt1)
  Ymax2 <- dnorm(Xmedia2, Xmedia2, Xdt2)

  ggplot(dat, aes(x)) +
    geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
    geom_ribbon(aes(ymin = y1, ymax = y2),
      data = subset(dat, x >= Xmedia2 + 1.5 * Xdt2),
      fill = "red", alpha = 0.5
    ) +
    annotate(
      geom = "segment",
      x = c(Xmedia1, Xmedia2), y = 0,
      xend = c(Xmedia1, Xmedia2), yend = c(Ymax1, Ymax2),
      linetype = "dashed",
      linewidth = 1,
      colour = c("grey", "black")
    ) +
    geom_line(aes(y = y1), linewidth = 1, colour = "grey") +
    geom_line(aes(y = y2), linewidth = 1, colour = "black") +
    theme(
      line = element_blank(),
      axis.line.y = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      axis.title = element_blank(),
      legend.position = "none",
      panel.grid = element_blank(),
      panel.background = element_rect(fill = "lightgray", colour = NA),
    )
}

graf_normal(250, 7, 253, 7)

Большое спасибо, Стефан, ваш код добавляет несколько хороших идей, таких как аннотация, лучше закодированная. Однако использование другой библиотеки не идеально для среды, в которой я должен развернуть решение. Я подумаю о вашем решении и r2evans.

— 19.02.2023 19:50

Привет Хуан. Конечно. Не стесняйтесь выбирать любой подход, который лучше всего соответствует вашим потребностям. И, конечно же, вы могли бы смешать оба, то есть вы могли бы заменить ggh4x::stat_difference в моем коде подходом geom_polygon от @r2evans.

— 19.02.2023 19:59

Я хотел отредактировать, но не могу через 5 минут, лол - ваш код добавляет МНОГО хороших идей!

— 19.02.2023 20:02

Стефан, я только что проверил, что ggh4x::stat_difference( можно напрямую заменить на geom_ribbon(, и код работает точно так же.

— 19.02.2023 23:40

Кстати, не могли бы вы объяснить использование тильды в ~ subset(.x, x >= Xmedia2 + 1.5 * Xdt2)? ничего не могу найти по этому поводу

— 20.02.2023 00:18

Ха-ха. Да, конечно. Думал слишком сложно. В вашем случае мы могли бы просто использовать geom_ribbon. stat_difference нужен только в том случае, если мы хотим заполнить, начиная с точки пересечения, или в случае, если мы хотим заполнить область между плотностями по всему диапазону данных. Относительно ~ subset(.x, x >= Xmedia2 + 1.5 * Xdt2): он подмножает данные, переданные ggplot(). В вашем случае мы могли бы просто сделать subset(dat, ...). Но использование этого обозначения удобно в случае, когда мы делаем что-то вроде dat %>% filter(...) %>% ggplot(...).

— 20.02.2023 01:42

Что ж, поверите или нет, вы предвосхитили мой следующий вопрос: ¿как заправиться с пересечения? Просто не делая подмножество и используя весь набор данных dat: вуаля, готово! Это весело :)

— 20.02.2023 10:02

Я сомневаюсь, что речь идет об обозначении тильды, с использованием тильды перед ~subset, а не с функцией subset - я никогда не использовал тильду перед функцией, и мне не ясно, для чего она используется. Спасибо!

— 20.02.2023 10:15

Привет Хуан. Вы найдете обозначение ~ повсюду в tidyverse. В основном ~ f(.x) то же самое, что и function(.x) f(.x), только короче. См., например. stackoverflow.com/questions/67643144/… некоторые ссылки.

— 20.02.2023 13:24

19.02.2023 19:40

Один из способов — использовать geom_polygon вместо stat_function.

Добавьте это в свою функцию перед ggplot():

  poly <- data.frame(xs = seq(Xmedia2+1.5*Xdt2, Xmax2, length.out = n)) |>
    transform(
      y1 = dnorm(xs, Xmedia1, Xdt1),
      y2 = dnorm(xs, Xmedia2, Xdt2)
    ) |>
    with(data.frame(X = c(xs, rev(xs)), Y = c(y1, rev(y2))))

А затем замените оба ваших вторых двух вызова stat_function одним

  geom_polygon(aes(X, Y), data = poly,
               fill = "red", alpha = 0.5) +

Полный источник:

graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2, n = 20) {

  Xmin1 <- Xmedia1-4*Xdt1
  Xmax1 <- Xmedia1+4*Xdt1

  Xmin2 <- Xmedia2-4*Xdt2
  Xmax2 <- Xmedia2+4*Xdt2

  Ymax1 <- max(dnorm(Xmedia1, Xmedia1, Xdt1))
  Ymax2 <- max(dnorm(Xmedia2, Xmedia2, Xdt2))

  Xmin <- min(Xmin1, Xmin2)
  Xmax <- max(Xmax1, Xmax2)

  poly <- data.frame(xs = seq(Xmedia2+1.5*Xdt2, Xmax2, length.out = n)) |>
    transform(
      y1 = dnorm(xs, Xmedia1, Xdt1),
      y2 = dnorm(xs, Xmedia2, Xdt2)
    ) |>
    with(data.frame(X = c(xs, rev(xs)), Y = c(y1, rev(y2))))

  ggplot(data.frame(X = c(Xmin, Xmax)), aes(x = X)) +
  geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
  geom_polygon(aes(X, Y), data = poly,
               fill = "red", alpha = 0.5) +
  stat_function(fun = dnorm, 
              args = c(Xmedia1, Xdt1), 
              linewidth = 1, 
              colour = "grey") +
  stat_function(fun = dnorm, 
              args = c(Xmedia2, Xdt2), 
              linewidth = 1, 
              colour = "black") +
  geom_segment(aes(x = Xmedia1, y = 0, xend = Xmedia1, yend = Ymax1), 
             linetype = "dashed", 
             linewidth = 0, 
             colour = "grey") +
  geom_segment(aes(x = Xmedia2, y = 0, xend = Xmedia2, yend = Ymax2), 
             linetype = "dashed", 
             linewidth = 0, 
             colour = "black") +
  ####################################################################
  # stat_function(fun = dnorm,
  #           args = c(Xmedia2, Xdt2),
  #           xlim = c(Xmedia2+1.5*Xdt2,Xmax2),
  #           geom = "area",
  #           fill = "red",
  #           alpha = 0.5) +
  # stat_function(fun = dnorm,
  #             args = c(Xmedia1, Xdt1),
  #             xlim = c(Xmedia1,Xmax1),
  #             geom = "area",
  #             fill = NA,
  #             alpha = 0.01) +
  ##################################################################
  theme(
    line = element_blank(),
    axis.line.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend.position = "none",
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "lightgray", colour = NA),
    ) +
    xlim(c(Xmin, Xmax)) 

}

Примечание. Я решил переместить geom_polygon намного раньше в стеке сюжетов, чтобы dnorm-строки были «поверх» красной области. Насколько это важно, зависит от вашего контекста и механизма рендеринга.

Мне нравится в этом решении, что код менее изменен, даже если мой код не слишком эффективен, я думаю, что он может быть понятнее в контексте, в котором я буду его использовать, когда люди учатся программировать. Из обоих ответов я вижу, что мне нужно сгенерировать значения Y, и я не могу полагаться на включенные функции. В любом случае не большая проблема. Большое спасибо, r2evans.

— 19.02.2023 19:53

Я не знаю ни одной функции, которая автоматизирует «область между двумя кривыми», кроме, возможно, предложения Стефана для ggh4x::stat_difference (я не пробовал). Рад, что это работает для вас!

— 19.02.2023 20:29

19.02.2023 19:40