В 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)

Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
3
0
50
2
Перейти к ответу Данный вопрос помечен как решенный

Ответы 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.

Juan Riera 19.02.2023 19:50

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

stefan 19.02.2023 19:59

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

Juan Riera 19.02.2023 20:02

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

Juan Riera 19.02.2023 23:40

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

Juan Riera 20.02.2023 00:18

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

stefan 20.02.2023 01:42

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

Juan Riera 20.02.2023 10:02

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

Juan Riera 20.02.2023 10:15

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

stefan 20.02.2023 13:24

Один из способов — использовать 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.

Juan Riera 19.02.2023 19:53

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

r2evans 19.02.2023 20:29

Другие вопросы по теме