{"spec_id":"probability-weibull","library":"ggplot2","language":"r","code":"#' anyplot.ai\n#' probability-weibull: Weibull Probability Plot for Reliability Analysis\n#' Library: ggplot2 3.5.1 | R 4.4.1\n#' Quality: 90/100 | Created: 2026-06-07\n\nlibrary(ggplot2)\nlibrary(scales)\nlibrary(ragg)\n\nset.seed(42)\n\n# Theme tokens\nTHEME       <- Sys.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG     <- if (THEME == \"light\") \"#FAF8F1\" else \"#1A1A17\"\nELEVATED_BG <- if (THEME == \"light\") \"#FFFDF6\" else \"#242420\"\nINK         <- if (THEME == \"light\") \"#1A1A17\" else \"#F0EFE8\"\nINK_SOFT    <- if (THEME == \"light\") \"#4A4A44\" else \"#B8B7B0\"\nINK_MUTED   <- if (THEME == \"light\") \"#6B6A63\" else \"#A8A79F\"\n\nIMPRINT_PALETTE <- c(\n  \"#009E73\",  # 1 — brand green (failures)\n  \"#AE3030\",  # 5 — matte red (censored — semantic: suspended/lost)\n  \"#4467A3\"   # 3 — blue (fitted line)\n)\n\n# Data: turbine blade fatigue-life data (hours to failure)\n# Simulate Weibull(shape=2.3, scale=8000) with some right-censored observations\nn_total <- 40\ntrue_shape <- 2.3\ntrue_scale <- 8000\n\n# Generate failure times\nall_times <- true_scale * (-log(runif(n_total)))^(1 / true_shape)\n# Right-censor at 10000 hours\ncensor_time <- 10000\nis_censored <- all_times > censor_time\ntime_to_failure <- pmin(all_times, censor_time)\n\n# Sort by time for proper rank ordering\nord <- order(time_to_failure)\ntime_to_failure <- time_to_failure[ord]\nis_censored <- is_censored[ord]\n\n# Count actual failures for median rank calculation\nn_fail <- sum(!is_censored)\nfailure_idx <- which(!is_censored)\n\n# Median rank (Benard's approximation): (i - 0.3) / (n + 0.4)\n# Only assign plotting positions to failures; censored are plotted separately\nmedian_ranks <- (seq_along(failure_idx) - 0.3) / (n_total + 0.4)\n\n# Weibull linearized y-axis: y = ln(-ln(1 - F))\nweibull_y <- log(-log(1 - median_ranks))\n\n# Data frame for failure points\ndf_fail <- data.frame(\n  time  = time_to_failure[failure_idx],\n  wb_y  = weibull_y,\n  type  = \"Failure\"\n)\n\n# For censored points: use the last assigned median rank as a ceiling (right-censored)\n# Place censored obs at a y that reflects \"at least survived this long\"\n# — just plot them at the axis mid-level annotated as censored; no formal rank\ncensored_times <- time_to_failure[is_censored]\ndf_cens <- data.frame(\n  time = censored_times,\n  wb_y = rep(median(weibull_y), length(censored_times)),\n  type = \"Censored (suspended)\"\n)\n\n# Fit line via OLS on linearised data (log(t) ~ wb_y)\nlog_times <- log(df_fail$time)\nfit <- lm(weibull_y ~ log_times, data = df_fail)\nbeta_hat  <- coef(fit)[2]            # slope = shape parameter\neta_hat   <- exp(-coef(fit)[1] / coef(fit)[2])  # scale parameter (characteristic life)\n\n# Fitted line over the data range\nx_seq     <- seq(log(min(df_fail$time) * 0.7), log(max(df_fail$time) * 1.3), length.out = 200)\ny_fit_seq <- coef(fit)[1] + coef(fit)[2] * x_seq\ndf_line   <- data.frame(time = exp(x_seq), wb_y = y_fit_seq)\n\n# 63.2% reference line: wb_y at F=0.632 = ln(-ln(1-0.632)) ≈ 0\nref_y  <- log(-log(1 - 0.632))  # ≈ -0.0006 ≈ 0\nref_63 <- data.frame(\n  x_start = min(df_line$time),\n  x_end   = exp(-coef(fit)[1] / coef(fit)[2]),  # = eta_hat\n  y_val   = ref_y\n)\n\n# Title\nplot_title <- \"probability-weibull · r · ggplot2 · anyplot.ai\"\ntitle_len  <- nchar(plot_title)\nbase_title_size <- 12\ntitle_size <- max(8, round(base_title_size * 67 / title_len))\n\n# Custom Weibull probability y-axis breaks and labels\n# Convert F values to wb_y = ln(-ln(1-F))\nf_breaks  <- c(0.01, 0.05, 0.10, 0.20, 0.30, 0.50, 0.632, 0.80, 0.90, 0.95, 0.99)\nwb_breaks <- log(-log(1 - f_breaks))\nf_labels  <- paste0(formatC(f_breaks * 100, format = \"g\", digits = 3), \"%\")\n\n# Plot\np <- ggplot() +\n  # Fitted Weibull line\n  geom_line(\n    data = df_line,\n    aes(x = time, y = wb_y),\n    color     = IMPRINT_PALETTE[3],\n    linewidth = 1.0,\n    linetype  = \"solid\"\n  ) +\n  # 63.2% horizontal reference line\n  geom_segment(\n    aes(x = ref_63$x_start, xend = ref_63$x_end,\n        y = ref_63$y_val,   yend = ref_63$y_val),\n    color     = INK_MUTED,\n    linewidth = 0.5,\n    linetype  = \"dashed\"\n  ) +\n  # Vertical drop to x-axis at eta\n  geom_segment(\n    aes(x = ref_63$x_end, xend = ref_63$x_end,\n        y = ref_63$y_val, yend = min(wb_breaks) - 0.3),\n    color     = INK_MUTED,\n    linewidth = 0.5,\n    linetype  = \"dashed\"\n  ) +\n  # Failure points (filled)\n  geom_point(\n    data = df_fail,\n    aes(x = time, y = wb_y, shape = \"Failure\", color = \"Failure\"),\n    size  = 2.8,\n    alpha = 0.85,\n    fill  = IMPRINT_PALETTE[1]\n  ) +\n  # Censored points (hollow)\n  geom_point(\n    data = df_cens,\n    aes(x = time, y = wb_y, shape = \"Censored (suspended)\", color = \"Censored (suspended)\"),\n    size   = 2.8,\n    alpha  = 0.85,\n    stroke = 1.2\n  ) +\n  # Parameter annotation\n  annotate(\n    \"label\",\n    x     = min(df_fail$time) * 1.1,\n    y     = max(wb_breaks) - 0.2,\n    label = sprintf(\"β (shape) = %.2f\\nη (scale) = %.0f h\", beta_hat, eta_hat),\n    hjust = 0,\n    vjust = 1,\n    size  = 3.5,\n    color = INK_SOFT,\n    fill  = ELEVATED_BG,\n    label.size = 0.3\n  ) +\n  # 63.2% label\n  annotate(\n    \"text\",\n    x     = ref_63$x_start * 1.05,\n    y     = ref_63$y_val + 0.08,\n    label = \"63.2%\",\n    hjust = 0,\n    size  = 2.5,\n    color = INK_MUTED\n  ) +\n  scale_x_log10(\n    labels = label_comma(),\n    name   = \"Time to Failure (hours)\"\n  ) +\n  scale_y_continuous(\n    breaks = wb_breaks,\n    labels = f_labels,\n    limits = c(min(wb_breaks) - 0.3, max(wb_breaks) + 0.1),\n    name   = \"Cumulative Failure Probability\"\n  ) +\n  scale_color_manual(\n    name   = NULL,\n    values = c(\"Failure\" = IMPRINT_PALETTE[1], \"Censored (suspended)\" = IMPRINT_PALETTE[2])\n  ) +\n  scale_shape_manual(\n    name   = NULL,\n    values = c(\"Failure\" = 19, \"Censored (suspended)\" = 1)\n  ) +\n  labs(title = plot_title) +\n  theme_minimal(base_size = 8) +\n  theme(\n    plot.background   = element_rect(fill = PAGE_BG, color = PAGE_BG),\n    panel.background  = element_rect(fill = PAGE_BG, color = NA),\n    panel.grid.major  = element_line(color = INK_SOFT, linewidth = 0.15),\n    panel.grid.minor  = element_blank(),\n    panel.border      = element_rect(color = INK_SOFT, fill = NA, linewidth = 0.4),\n    axis.title        = element_text(color = INK,      size = 10),\n    axis.text         = element_text(color = INK_SOFT, size = 8),\n    axis.text.x       = element_text(color = INK_SOFT, size = 8),\n    plot.title        = element_text(color = INK,      size = title_size, face = \"bold\",\n                                     margin = margin(b = 8)),\n    legend.background = element_rect(fill = ELEVATED_BG, color = INK_SOFT, linewidth = 0.3),\n    legend.text       = element_text(color = INK_SOFT, size = 8),\n    legend.title      = element_text(color = INK,      size = 9),\n    legend.position   = \"bottom\",\n    legend.key        = element_rect(fill = NA, color = NA),\n    plot.margin       = margin(12, 16, 12, 12)\n  )\n\n# Save\nggsave(\n  filename = sprintf(\"plot-%s.png\", THEME),\n  plot     = p,\n  device   = ragg::agg_png,\n  width    = 8,\n  height   = 4.5,\n  units    = \"in\",\n  dpi      = 400\n)\n"}