Please update this code to have 4 cutoff points and 5 ratings:


# Ordinal intuition
library(tidyverse)

d <- 
  tibble(z = seq(from = -3, to = 3, by = .1)) %>% 
  mutate(`p(z)`   = dnorm(z, mean = 0, sd = 1),
         `Phi(z)` = pnorm(z, mean = 0, sd = 1))

library(scico)

sl <- scico(palette = "lajolla", n = 9)

scales::show_col(sl)

theme_set(
  theme_linedraw() +
    theme(panel.grid = element_blank(),
          strip.background = element_rect(color = sl[9], fill = sl[9]),
          strip.text = element_text(color = sl[1]))
)

p1 <-
  d %>% 
  ggplot(aes(x = z, y = `p(z)`)) +
  geom_area(aes(fill = z <= 1),
            show.legend = F) +
  geom_line(linewidth = 1, color = sl[3]) +
  scale_fill_manual(values = c("transparent", sl[2])) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + 
  labs(title = "Normal Density",
       y = expression(p(italic(z))))

p2 <-
  d %>% 
  ggplot(aes(x = z, y = `Phi(z)`)) +
  geom_area(aes(fill = z <= 1),
            show.legend = F) +
  geom_line(linewidth = 1, color = sl[3]) +
  scale_fill_manual(values = c("transparent", sl[2])) +
  scale_y_continuous(expand = expansion(mult = c(0, 0)), limits = 0:1) + 
  labs(title = "Cumulative Normal",
       y = expression(Phi(italic(z))))

library(patchwork)

p1 / p2 & 
  scale_x_continuous(breaks = -2:2) & 
  coord_cartesian(xlim = c(-2.5, 2.5))

(theta_1 <- seq(from = 1.5, to = 6.5, by = 1))
(label_1 <- 1:7)

den <-
  # define the parameters for the underlying normal distribution
  tibble(mu    = 4,
         sigma = 1.5) %>% 
  mutate(strip = str_c("mu==", mu, "~~sigma==", sigma)) %>% 
  # this will allow us to rescale the density in terms of the bar plot
  mutate(multiplier = 26 / dnorm(mu, mu, sigma)) %>% 
  # we need values for the x-axis
  expand_grid(y = seq(from = -1, to = 9, by = .1)) %>% 
  # compute the density values
  mutate(density = dnorm(y, mu, sigma)) %>% 
  # use that multiplier value from above to rescale the density values
  mutate(percent = density * multiplier)

head(den)

bar <-
  # define the parameters for the underlying normal distribution
  tibble(mu    = 4,
         sigma = 1.5) %>% 
  mutate(strip = str_c("mu==", mu, "~~sigma==", sigma)) %>% 
  # take random draws from the underlying normal distribution
  mutate(draw = map2(mu, sigma, ~rnorm(1e4, mean = .x, sd = .y))) %>% 
  unnest(draw) %>% 
  # bin those draws into ordinal categories defined by `theta_1`
  # and named by `label_1`
  mutate(y = case_when(
    draw  < theta_1[1] ~ label_1[1],
    draw  < theta_1[2] ~ label_1[2],
    draw  < theta_1[3] ~ label_1[3],
    draw  < theta_1[4] ~ label_1[4],
    draw  < theta_1[5] ~ label_1[5],
    draw  < theta_1[6] ~ label_1[6],
    draw >= theta_1[6] ~ label_1[7]
  )) %>% 
  # summarize
  count(y) %>% 
  mutate(percent = (100 * n / sum(n)) %>% round(0))   %>% 
  mutate(percent_label = str_c(percent, "%"),
         percent_max   = max(percent))

head(bar)

bar %>% 
  ggplot(aes(x = y, y = percent)) +
  geom_area(data = den,
            fill = sl[3]) +
  geom_vline(xintercept = theta_1, color = "white", linetype = 3) +
  geom_col(width = .5, alpha = .85, fill = sl[7]) +
  geom_text(aes(y = percent + 2, label = percent_label),
            size = 3.5) +
  annotate(geom = "text", 
           x = theta_1, y = -5.55, label = theta_1, 
           size = 3) +
  scale_x_continuous(NULL, expand = c(0, 0),
                     breaks = theta_1,
                     labels = parse(text = str_c("theta[", 1:6, "]"))) +
  scale_y_continuous(NULL, breaks = NULL, expand = expansion(mult = c(0, 0.06))) +
  coord_cartesian(ylim = c(0, 28.5), clip = F) +
  theme(plot.margin = margin(5.5, 5.5, 11, 5.5)) +
  facet_wrap(~ strip, labeller = label_parsed)