How to plot multiple functions with stat_function()

Prep steps

penguins <- palmerpenguins::penguins |> 
  filter(!is.na(sex))

grouped_penguins <-  penguins |> 
  mutate(species = if_else(species == 'Gentoo', species, 'Others'))

penguin_colors <- c("#E69F00", "#009E73", "#0072B2", 'grey80')
names(penguin_colors) <- c(
  unique(penguins$species |> as.character()), 
  'Others'
)

density_chart <- grouped_penguins |> 
  ggplot() +
  geom_density(
    aes(x = body_mass_g, fill = species), 
    alpha = 0.8,
    col = NA
  ) +
  labs(
    x = 'Body weight (in g)',
    y = element_blank(),
    title = 'Penguin weights are somewhat normally distributed',
    fill = 'Species'
  ) +
  scale_fill_manual(values = penguin_colors) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.025)))

params_grouped_penguins <- grouped_penguins |> 
  summarize(
    mean_body_mass_g = mean(body_mass_g),
    sd_body_mass_g   = sd(body_mass_g),
    .by              = species
  )

Add function

# Get parameters
params_grouped_penguins
## # A tibble: 2 × 3
##   species mean_body_mass_g sd_body_mass_g
##   <chr>              <dbl>          <dbl>
## 1 Others             3715.           436.
## 2 Gentoo             5092.           501.


# Compute list of layers with pmap()
fun_layers <- params_grouped_penguins |> 
  pmap(
    ~{
      # Pass current parameters to new list
        l <- list(...)
        
      # Use list to put parameters into correct place in dnorm()
        stat_function(
          fun = \(x) dnorm(
            x, 
            mean = l$mean_body_mass_g,
            sd = l$sd_body_mass_g
          ),
          
          # Add some aesthetics as well
          linewidth = 1,
          linetype = if (l$species == 'Gentoo') 2 else 3
        )
      }
  )

# Add layers to a previously computed ggplot
density_chart +
  fun_layers