Before Seeing Data

p_1 <- tibble(x_canvas = c(0, 1))

alpha <- 1
beta <- 1

plot_p1_fcn <- function(a, b) {
  p_1 %>%
    ggplot(aes(x = x_canvas)) +
    stat_function(
      fun = dbeta,
      args = list(a, b),
      color = "#2c3e50",
      size = 1,
      alpha = .8
    ) +
    #  ylim(c(0, 1.5)) +
    labs(
      y = "",
      x = "P1: Probability of Fracture",
      title = "Credible Estimates for P1",
      subtitle = str_glue("Beta ({a}, {b})")
    ) +
    theme_bw() +
    theme(
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank()
    )
}

plot_p1_fcn(a = alpha, b = beta)

Add Round 1 Benchtop Data

Recall from above:

Round 1: 6 parts tested and all 6 passed (0 failures). So we add 6 to the cumulative sum for b and nothing to a.

alpha <- 1
beta <- 7

plot_p1_fcn(a = alpha, b = beta)

Add Round 2 Benchtop Data

Recall from above:

Round 2: 5 parts tested and all 5 passed (0 failures). So we add 5 to the cumulative sum for b and nothing to a.

alpha <- 1
beta <- 12

plot_p1_fcn(a = alpha, b = beta)

Add Round 3 Benchtop Data

Recall from above:

Round 3: 12 parts tested and all 11 passed (1 failure). So we add 11 to the cumulative sum for b and 1 to the cumulative sum for a.

alpha <- 2
beta <- 23

plot_p1_fcn(a = alpha, b = beta)

Pausing here to assess, we can see a wide range of values for P1 that are still consistent with the data, ranging from close to 0 up to around .25.

Prior Probability based on FEA (Input 1)

p_2 <- tibble(x_canvas = c(0, 1))

alpha <- 1
beta <- 20

plot_p2_fcn <- function(a, b) {
  p_2 %>%
    ggplot(aes(x = x_canvas)) +
    stat_function(
      fun = dbeta,
      args = list(a, b),
      color = "firebrick",
      size = 1,
      alpha = .8
    ) +
    #  ylim(c(0, 1.5)) +
    labs(
      y = "",
      x = "P2: Probability of Fracture resulting in Harm",
      title = "Credible Estimates for P2",
      subtitle = str_glue("Beta ({a}, {b})")
    ) +
    theme_bw() +
    theme(
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank()
    )
}

plot_p2_fcn(a = alpha, b = beta) +
  labs(caption = "prior probabilities assessed primarily from FEA")

Add Input 2 (Benchtop Data)

Now we add our study where we looked at 6 parts that were intentionally fractured and saw no propagation or cascading failures. So we add 6 to the cumulative sum for b and nothing to a.

alpha <- 1
beta <- 26

plot_p2_fcn(alpha, beta)

Sample from P1 Posterior and Visualize

set.seed(123)
n_draws <- 500000

p1_post_draws_tbl <- tibble(draws = rbeta(n = n_draws, shape1 = 2, shape2 = 23))

p1_post_draws_tbl %>%
  gt_preview()

p1_post_draws_tbl %>%
  ggplot(aes(x = draws)) +
  geom_histogram(
    binwidth = .002,
    color = "white",
    fill = "#2c3e50",
    alpha = 0.8
  ) +
  xlim(c(-0.05, .5)) +
  labs(
    y = "",
    title = "Samples from P1 Posterior",
    subtitle = "Current State of Belief",
    x = "P1: Probability of Fracture"
  ) +
  theme_bw() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

Sample from P2 Posterior and Visualize

set.seed(123)
n_draws <- 500000

p2_post_draws_tbl <- tibble(draws = rbeta(n = n_draws, shape1 = 1, shape2 = 26))

p2_post_draws_tbl %>%
  gt_preview()

p2_post_draws_tbl %>%
  ggplot(aes(x = draws)) +
  geom_histogram(
    binwidth = .002,
    color = "white",
    fill = "firebrick",
    alpha = 0.8
  ) +
  xlim(c(-0.05, .5)) +
  labs(
    y = "",
    title = "Samples from P2 Posterior",
    subtitle = "Current State of Belief",
    x = "P2: Probability of a Fracture Resulting in a Harm"
  ) +
  theme_bw() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

Combine to Estimate Ph

joint_post_tbl <- p1_post_draws_tbl %>%
  bind_cols(p2_post_draws_tbl) %>%
  rename(
    p1_draws = draws...1,
    p2_draws = draws...2
  ) %>%
  mutate(p1_x_p2 = p1_draws * p2_draws)

joint_post_tbl %>%
  ggplot(aes(x = p1_x_p2)) +
  geom_histogram(
    binwidth = .0001,
    color = "white",
    fill = "purple",
    alpha = 0.8
  ) +
  labs(
    y = "",
    title = "Ph - Probabilty Model of Fracture Occurring AND Resulting in Harm",
    subtitle = "Current State",
    x = "Ph: Probabilty of Fracture Occurring AND Resulting in Harm"
  ) +
  theme_bw() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  xlim(c(-.001, .03))

The above plot can also be summarized in tabular form.

joint_post_tbl %>%
  select(p1_x_p2) %>%
  mutate_if(is.numeric, round, 4) %>%
  median_qi(.width = .95) %>%
  gt() %>%
  cols_align("center") %>%
  opt_row_striping(row_striping = TRUE) %>%
  tab_header(
    title = "Summary of Ph Posterior",
    subtitle = "Based on Bayesian Beta/Binomial Model"
  ) %>%
  cols_width(everything() ~ px(100)) %>%
  tab_style(
    style = list(
      cell_text(weight = "bold")
    ),
    locations = cells_column_labels(columns = everything())
  )

We have arrived at an estimate for Ph, the probability of the fracture occurring and resulting in a harm. A point estimate is .0016 or 16 out of 1000. A decision made off this single value may result in pushing a product that seems relatively safe. But the 95% credible interval covers up to .014 or 1.4 in 100 which would not be acceptable in most circumstances. The Bayesian approach has allowed us to carry that uncertainty all the way through the calculations to make a more informed decision.

This is the distillation of a lot of work to prepare the problem for a business decision: do nothing, initiate a redesign, or proceed with more testing.

In this hypothetical example, credible risk values above 1% are likely too high to proceed without design changes or additional testing.

If you’ve made it this far, I thank you for your attention.