|
| 1 | +""" pyplots.ai |
| 2 | +probability-weibull: Weibull Probability Plot for Reliability Analysis |
| 3 | +Library: pygal 3.1.0 | Python 3.14.3 |
| 4 | +Quality: 87/100 | Created: 2026-03-11 |
| 5 | +""" |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +import pygal |
| 9 | +from pygal.style import Style |
| 10 | +from scipy import stats |
| 11 | + |
| 12 | + |
| 13 | +# Data — turbine blade fatigue-life (hours) with failures and suspensions |
| 14 | +np.random.seed(42) |
| 15 | +n_failures = 18 |
| 16 | +n_censored = 5 |
| 17 | +beta_true = 2.5 |
| 18 | +eta_true = 8000 |
| 19 | + |
| 20 | +failure_times = np.sort(stats.weibull_min.rvs(beta_true, scale=eta_true, size=n_failures)) |
| 21 | +censored_times = np.sort(np.random.uniform(2000, 9000, n_censored)) |
| 22 | + |
| 23 | +all_times = np.concatenate([failure_times, censored_times]) |
| 24 | +is_failure = np.concatenate([np.ones(n_failures), np.zeros(n_censored)]) |
| 25 | + |
| 26 | +sort_idx = np.argsort(all_times) |
| 27 | +all_times = all_times[sort_idx] |
| 28 | +is_failure = is_failure[sort_idx] |
| 29 | + |
| 30 | +# Median rank plotting positions (i-0.3)/(n+0.4) for failures only |
| 31 | +failure_ranks = np.cumsum(is_failure) |
| 32 | +total_failures = failure_ranks[-1] |
| 33 | +median_ranks = (failure_ranks - 0.3) / (total_failures + 0.4) |
| 34 | + |
| 35 | +failure_mask = is_failure.astype(bool) |
| 36 | +failure_x = all_times[failure_mask] |
| 37 | +failure_prob = median_ranks[failure_mask] |
| 38 | + |
| 39 | +# Weibull linearization: ln(-ln(1-F)) for y-axis, ln(time) for x-axis |
| 40 | +weibull_y_failures = np.log(-np.log(1.0 - failure_prob)) |
| 41 | +ln_x_failures = np.log(failure_x) |
| 42 | + |
| 43 | +# Fit line using least squares on linearized data |
| 44 | +slope, intercept, r_value, _, _ = stats.linregress(ln_x_failures, weibull_y_failures) |
| 45 | +beta_est = slope |
| 46 | +eta_est = np.exp(-intercept / beta_est) |
| 47 | + |
| 48 | +# Fitted line spanning full data range |
| 49 | +x_fit_range = np.linspace(np.log(min(all_times) * 0.7), np.log(max(all_times) * 1.3), 100) |
| 50 | +y_fit_line = slope * x_fit_range + intercept |
| 51 | + |
| 52 | +# Censored points — place on the fitted line at their log-time |
| 53 | +censored_x = all_times[~failure_mask] |
| 54 | +ln_censored_x = np.log(censored_x) |
| 55 | +censored_y_on_line = slope * ln_censored_x + intercept |
| 56 | + |
| 57 | +# 63.2% reference line (characteristic life): F = 0.632 → ln(-ln(1-0.632)) ≈ 0.0 |
| 58 | +ref_y = np.log(-np.log(1 - 0.632)) |
| 59 | + |
| 60 | +# B10 life: F = 0.10 → time where 10% of units have failed |
| 61 | +b10_y = np.log(-np.log(1 - 0.10)) |
| 62 | +b10_ln_x = (b10_y - intercept) / slope |
| 63 | +b10_hours = np.exp(b10_ln_x) |
| 64 | + |
| 65 | +# Characteristic life intersection: where fitted line crosses 63.2% |
| 66 | +eta_ln_x = (ref_y - intercept) / slope |
| 67 | + |
| 68 | +# Axis labels — convert back from log/Weibull scale for readability |
| 69 | +x_tick_values = [1000, 2000, 3000, 5000, 7000, 10000, 15000] |
| 70 | +x_tick_ln = [np.log(v) for v in x_tick_values] |
| 71 | + |
| 72 | +prob_levels = [0.01, 0.05, 0.10, 0.20, 0.50, 0.632, 0.80, 0.90, 0.95, 0.99] |
| 73 | +y_tick_weibull = [np.log(-np.log(1.0 - p)) for p in prob_levels] |
| 74 | +y_tick_labels = [f"{p * 100:.1f}%" if p == 0.632 else f"{p * 100:.0f}%" for p in prob_levels] |
| 75 | + |
| 76 | +# Style — cohesive blue-tonal palette with warm accent for annotations |
| 77 | +font = "DejaVu Sans, Helvetica, Arial, sans-serif" |
| 78 | +custom_style = Style( |
| 79 | + background="white", |
| 80 | + plot_background="#fafbfc", |
| 81 | + foreground="#2c3e50", |
| 82 | + foreground_strong="#1a252f", |
| 83 | + foreground_subtle="#e0e4e8", |
| 84 | + guide_stroke_color="#e0e4e8", |
| 85 | + guide_stroke_dasharray="4, 8", |
| 86 | + colors=( |
| 87 | + "#4a7fb5", # fit line — medium blue |
| 88 | + "#1b3a5c", # failures — dark navy |
| 89 | + "#7cadd4", # censored — light steel blue |
| 90 | + "#95a5a6", # 63.2% reference — cool gray |
| 91 | + "#d35400", # η marker — warm burnt orange (accent) |
| 92 | + "#27ae60", # B10 marker — muted green (distinct from blue) |
| 93 | + ), |
| 94 | + font_family=font, |
| 95 | + title_font_family=font, |
| 96 | + title_font_size=56, |
| 97 | + label_font_size=40, |
| 98 | + major_label_font_size=36, |
| 99 | + legend_font_size=38, |
| 100 | + legend_font_family=font, |
| 101 | + value_font_size=32, |
| 102 | + tooltip_font_size=28, |
| 103 | + tooltip_font_family=font, |
| 104 | + opacity=0.92, |
| 105 | + opacity_hover=1.0, |
| 106 | + stroke_opacity=1, |
| 107 | + stroke_opacity_hover=1, |
| 108 | +) |
| 109 | + |
| 110 | +# Axis bounds — tighter y-range to reduce wasted space at extremes |
| 111 | +x_min_ln = np.log(800) |
| 112 | +x_max_ln = np.log(18000) |
| 113 | +y_min_w = np.log(-np.log(1 - 0.008)) |
| 114 | +y_max_w = np.log(-np.log(1 - 0.993)) |
| 115 | + |
| 116 | +# Chart configuration |
| 117 | +chart = pygal.XY( |
| 118 | + width=4800, |
| 119 | + height=2700, |
| 120 | + style=custom_style, |
| 121 | + title="Turbine Blade Fatigue Life · probability-weibull · pygal · pyplots.ai", |
| 122 | + x_title="Time to Failure (hours, log scale)", |
| 123 | + y_title="Cumulative Failure Probability", |
| 124 | + show_legend=True, |
| 125 | + legend_at_bottom=True, |
| 126 | + legend_at_bottom_columns=3, |
| 127 | + legend_box_size=30, |
| 128 | + stroke=False, |
| 129 | + dots_size=10, |
| 130 | + show_x_guides=True, |
| 131 | + show_y_guides=True, |
| 132 | + margin_bottom=100, |
| 133 | + margin_left=90, |
| 134 | + margin_right=60, |
| 135 | + margin_top=60, |
| 136 | + truncate_legend=-1, |
| 137 | + range=(y_min_w, y_max_w), |
| 138 | + xrange=(x_min_ln, x_max_ln), |
| 139 | + print_values=False, |
| 140 | + print_zeroes=False, |
| 141 | + js=[], |
| 142 | + x_labels=[float(v) for v in x_tick_ln], |
| 143 | + y_labels=[float(v) for v in y_tick_weibull], |
| 144 | + x_value_formatter=lambda x: f"{np.exp(x):,.0f}h", |
| 145 | + value_formatter=lambda y: f"{(1 - np.exp(-np.exp(y))) * 100:.1f}%", |
| 146 | + tooltip_border_radius=10, |
| 147 | + tooltip_fancy_mode=True, |
| 148 | + dynamic_print_values=True, |
| 149 | +) |
| 150 | + |
| 151 | +# Override label display for axes |
| 152 | +chart.y_labels = [{"value": float(v), "label": lbl} for v, lbl in zip(y_tick_weibull, y_tick_labels, strict=True)] |
| 153 | +chart.x_labels = [{"value": float(v), "label": f"{int(t):,}"} for v, t in zip(x_tick_ln, x_tick_values, strict=True)] |
| 154 | + |
| 155 | +# Fitted line — medium blue, prominent stroke |
| 156 | +fit_points = [(float(x), float(y)) for x, y in zip(x_fit_range, y_fit_line, strict=True)] |
| 157 | +chart.add( |
| 158 | + f"Weibull Fit (β={beta_est:.2f}, η={eta_est:,.0f}h, R²={r_value**2:.3f})", |
| 159 | + fit_points, |
| 160 | + stroke=True, |
| 161 | + show_dots=False, |
| 162 | + stroke_style={"width": 8, "linecap": "round", "linejoin": "round"}, |
| 163 | +) |
| 164 | + |
| 165 | +# Failure data points — dark navy with per-point tooltip labels |
| 166 | +failure_points = [ |
| 167 | + { |
| 168 | + "value": (float(x), float(y)), |
| 169 | + "label": f"Failure at {np.exp(x):,.0f}h — F={((1 - np.exp(-np.exp(y))) * 100):.1f}%", |
| 170 | + } |
| 171 | + for x, y in zip(ln_x_failures, weibull_y_failures, strict=True) |
| 172 | +] |
| 173 | +chart.add(f"Failures (n={n_failures})", failure_points, stroke=False, dots_size=12) |
| 174 | + |
| 175 | +# Censored data points — light steel blue, visually lighter than failures |
| 176 | +censored_points = [ |
| 177 | + {"value": (float(x), float(y)), "label": f"Censored at {np.exp(x):,.0f}h (suspended test)"} |
| 178 | + for x, y in zip(ln_censored_x, censored_y_on_line, strict=True) |
| 179 | +] |
| 180 | +chart.add(f"Censored (n={n_censored})", censored_points, stroke=False, dots_size=9) |
| 181 | + |
| 182 | +# 63.2% reference line — cool gray dashed |
| 183 | +ref_line = [(float(x_min_ln), float(ref_y)), (float(x_max_ln), float(ref_y))] |
| 184 | +chart.add( |
| 185 | + "63.2% Characteristic Life", |
| 186 | + ref_line, |
| 187 | + stroke=True, |
| 188 | + show_dots=False, |
| 189 | + stroke_style={"width": 4, "dasharray": "16, 10", "linecap": "round"}, |
| 190 | +) |
| 191 | + |
| 192 | +# η marker — warm burnt orange accent, on-chart printed value |
| 193 | +eta_marker = [ |
| 194 | + {"value": (float(eta_ln_x), float(ref_y)), "label": f"η = {eta_est:,.0f}h (Characteristic Life at 63.2%)"} |
| 195 | +] |
| 196 | +chart.add( |
| 197 | + f"η = {eta_est:,.0f}h", |
| 198 | + eta_marker, |
| 199 | + stroke=False, |
| 200 | + dots_size=22, |
| 201 | + print_values=True, |
| 202 | + formatter=lambda x: f"η = {eta_est:,.0f}h", |
| 203 | +) |
| 204 | + |
| 205 | +# B10 marker — muted green, distinct from blues, on-chart printed value |
| 206 | +b10_marker = [{"value": (float(b10_ln_x), float(b10_y)), "label": f"B10 = {b10_hours:,.0f}h (10% Failure Life)"}] |
| 207 | +chart.add( |
| 208 | + f"B10 = {b10_hours:,.0f}h", |
| 209 | + b10_marker, |
| 210 | + stroke=False, |
| 211 | + dots_size=22, |
| 212 | + print_values=True, |
| 213 | + formatter=lambda x: f"B10 = {b10_hours:,.0f}h", |
| 214 | +) |
| 215 | + |
| 216 | +# Save |
| 217 | +chart.render_to_png("plot.png") |
| 218 | +chart.render_to_file("plot.html") |
0 commit comments