Hierarchical Bayesian Analysis to Lead Decision-Making by Varying Factory and Line Variation

Title & Overview

Even factories or lines producing the same product can vary in defect rates, yields, and temperature sensitivity. In this article, focusing on a fictional precision parts manufacturer, we will focus on “Partial Pooling” Sharing Overall Trends While Leaving Differences Among Each Group and implement No.071 to No.080 in a single flow.

The question addressed is: “Can a small-scale production line with a high observed defect rate be immediately identified as a problem line?” Using a hierarchical structure, the strength of the estimate varies according to data volume, allowing for improvement, additional measurement, and horizontal deployment.

[!NOTE] This material is a notebook previously used by Surikoubo (or personally by the representative, Kazuyama), and has been reconstructed, edited, and published with the company’s permission. All data listed is fictional and has no relation whatsoever to real companies, factories, or figures.

Introduction: Practical Challenges in Manufacturing Covered in This Article

In multi-site, multi-product quality control, ranking the main lines with a large base and prototyping/small-batch lines on the same standard can make a chance shift look like a “difference in ability.” In this article, we estimate each line by borrowing information shared by the group and consider where to allocate the limited improvement resources.

Common situations on site

• Sites are evaluated based solely on monthly defect rate rankings.
• The newly installed lines have fewer inspections, with extreme values such as 0% or 10%.
• Varieties, materials, and equipment conditions are mixed, making simple averages impossible to compare
• There are two choices: “replace with company-wide average” or “treat each site completely differently.”

Why is this issue so difficult to judge?

Observed errors include both true process differences and sampling errors simultaneously. In particular, groups with fewer tests have larger sampling errors, and point estimation alone cannot distinguish them. Furthermore, ignoring field differences dilutes anomalies, and completely separating sites makes estimation unstable.

Overview of Exercise covered this time

No.071–080 cover hierarchical structure, three types of pooling, rate estimation, minority stabilization, random intercepts and coefficients, visualization, and introduction decision, in that order. The central beta and binomial model is

y_j\mid p_j\sim\mathrm{Binomial}(n_j,p_j),\qquad p_j\sim\mathrm{Beta}(lpha,eta)

That’s right. The posterior distribution is p_j\mid y_j\sim\mathrm{Beta}(lpha+y_j,eta+n_j-y_j), and the posterior mean is

E[p_j\mid y_j]=rac{n_j}{n_j+lpha+eta}rac{y_j}{n_j}+rac{lpha+eta}{n_j+lpha+eta}rac{lpha}{lpha+eta}

That’s right. In other words, groups with a lot of data tend to place more emphasis on actual results, while groups with less data tend to emphasize overall trends.

Preparing the Python environment

No external data is used; random number seeds are fixed. scipy is used for distribution and interval estimation, and matplotlib is used for visualization.

import sys
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from scipy import stats

SEED = 20260712
rng = np.random.default_rng(SEED)
pd.set_option("display.precision", 4)
print(f"Python {sys.version.split()[0]}")
print(f"NumPy {np.__version__} / pandas {pd.__version__} / matplotlib {matplotlib.__version__}")

Python 3.13.1
NumPy 2.5.1 / pandas 3.0.3 / matplotlib 3.11.0

Creation of Fictional Data

8. Generates inspection and defect numbers for production lines, 6. Number of product passes, improvement responses by material segment, and daily process data by line. The true value is set to show the difference, but true values are not used for analysis.

lines = [f"Line-{c}" for c in "ABCDEFGH"]
n_inspected = np.array([1200, 900, 650, 400, 180, 80, 35, 12])
true_defect = np.array([0.018, 0.022, 0.027, 0.020, 0.035, 0.025, 0.045, 0.030])
defects = rng.binomial(n_inspected, true_defect)
line_df = pd.DataFrame({"line": lines, "inspected": n_inspected, "defects": defects})
line_df["observed_rate"] = line_df["defects"] / line_df["inspected"]

products = [f"Product-{i}" for i in range(1, 7)]
prod_n = np.array([1000, 700, 350, 160, 60, 20])
true_yield = np.array([0.965, 0.952, 0.970, 0.940, 0.955, 0.930])
passed = rng.binomial(prod_n, true_yield)
product_df = pd.DataFrame({"product": products, "produced": prod_n, "passed": passed})
product_df["observed_yield"] = product_df["passed"] / product_df["produced"]

line_df

	line	inspected	defects	observed_rate
0	Line-A	1200	24	0.0200
1	Line-B	900	19	0.0211
2	Line-C	650	22	0.0338
3	Line-D	400	7	0.0175
4	Line-E	180	2	0.0111
5	Line-F	80	1	0.0125
6	Line-G	35	4	0.1143
7	Line-H	12	1	0.0833

No.071: Understanding the Hierarchical Bayesian Concept

Meaning in Practice

Although the lines operate independently, they share the same equipment standards, material specifications, and quality systems. This “different but not irrelevant” expression is the hierarchical model.

Approach to Analysis and Modeling

Here, the common population is $\mathrm{Beta}(1,39)$, with an average company-wide defect rate of 2.5% and a pre-information volume equivalent to 40 cases. This is not a value to fix a conclusion, but rather an assumption that should be analyzed for sensitivity later.

Check with Python

alpha0, beta0 = 1.0, 39.0
line_df["post_alpha"] = alpha0 + line_df["defects"]
line_df["post_beta"] = beta0 + line_df["inspected"] - line_df["defects"]
line_df["hier_rate"] = line_df["post_alpha"] / (line_df["post_alpha"] + line_df["post_beta"])
line_df[["line", "inspected", "defects", "observed_rate", "hier_rate"]]

	line	inspected	defects	observed_rate	hier_rate
0	Line-A	1200	24	0.0200	0.0202
1	Line-B	900	19	0.0211	0.0213
2	Line-C	650	22	0.0338	0.0333
3	Line-D	400	7	0.0175	0.0182
4	Line-E	180	2	0.0111	0.0136
5	Line-F	80	1	0.0125	0.0167
6	Line-G	35	4	0.1143	0.0667
7	Line-H	12	1	0.0833	0.0385

Reading the results

In lines with a high number of tests, the estimated hierarchy values are close to the actual rate, while in the small groups, they strongly approach 2.5%. This move is not a correction to hide the data, but rather a reflection of the magnitude of the sampling error in the estimate.

No.072: Comparing Complete, Non-Pooling, and Partial Pooling

Meaning in Practice

Full pooling treats all lines as the same, while non-pooling treats each line as if it were a separate company. Partial pooling shares information according to the amount of data in between.

Approach to Analysis and Modeling

Complete pooling refers to the overall defect rate, non-pooling is the observed rate, and partial pooling is the post-mortem average mentioned earlier. Listing these three options explains the assumptions implicitly set by the evaluation system.

Check with Python

pooled_rate = line_df["defects"].sum() / line_df["inspected"].sum()
comparison = line_df[["line", "inspected", "observed_rate", "hier_rate"]].copy()
comparison["complete_pooling"] = pooled_rate
comparison = comparison.rename(columns={"observed_rate": "no_pooling", "hier_rate": "partial_pooling"})
comparison[["line", "inspected", "complete_pooling", "no_pooling", "partial_pooling"]]

	line	inspected	complete_pooling	no_pooling	partial_pooling
0	Line-A	1200	0.0231	0.0200	0.0202
1	Line-B	900	0.0231	0.0211	0.0213
2	Line-C	650	0.0231	0.0338	0.0333
3	Line-D	400	0.0231	0.0175	0.0182
4	Line-E	180	0.0231	0.0111	0.0136
5	Line-F	80	0.0231	0.0125	0.0167
6	Line-G	35	0.0231	0.1143	0.0667
7	Line-H	12	0.0231	0.0833	0.0385

plot_df = comparison.set_index("line")[["complete_pooling", "no_pooling", "partial_pooling"]]
ax = plot_df.plot(kind="bar", figsize=(10, 4))
ax.set_title("Three pooling strategies for defect rates")
ax.set_xlabel("Production line")
ax.set_ylabel("Defect rate")
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

Reading the results

In full pooling, differences in areas for improvement disappear, while in non-pooling, extremes of a minority group are emphasized. Partial pooling is suitable for primary screening because it suppresses extreme ranking fluctuations while leaving a difference.

No.073: Estimating Store-Specific CVR Using Hierarchical Bayes

Meaning in Practice

The store-specific CVR in the original title is binary data for manufacturing as ‘defect occurrence rate per inspection opportunity.’ It can be converted into a quality rate by site or line.

Approach to Analysis and Modeling

From the posterior distribution of each line, we calculate the 95% credit range and the “probability that the non-performing ratio exceeds 3%.” Priority is determined not by point estimation, but by the probability of exceeding quality standards.

Check with Python

line_df["ci_low"] = stats.beta.ppf(0.025, line_df["post_alpha"], line_df["post_beta"])
line_df["ci_high"] = stats.beta.ppf(0.975, line_df["post_alpha"], line_df["post_beta"])
line_df["prob_over_3pct"] = 1 - stats.beta.cdf(0.03, line_df["post_alpha"], line_df["post_beta"])
line_df[["line", "hier_rate", "ci_low", "ci_high", "prob_over_3pct"]].sort_values("prob_over_3pct", ascending=False)

	line	hier_rate	ci_low	ci_high	prob_over_3pct
6	Line-G	0.0667	0.0223	0.1327	0.9283
2	Line-C	0.0333	0.0213	0.0479	0.6692
7	Line-H	0.0385	0.0048	0.1045	0.5452
5	Line-F	0.0167	0.0020	0.0459	0.1248
3	Line-D	0.0182	0.0079	0.0326	0.0470
1	Line-B	0.0213	0.0131	0.0314	0.0425
4	Line-E	0.0136	0.0028	0.0326	0.0388
0	Line-A	0.0202	0.0131	0.0287	0.0133

Reading the results

A high prob_over_3pct line is not just a high observation rate, but a candidate with a high likelihood of exceeding the standard. For lines with a wide credit range, additional inspections are prioritized over definitive assessments.

No.074: Estimating Sales Rates by Product Using Hierarchical Bayes

Meaning in Practice

The original title translates the sales rate into product-specific yield. Even when the number of parameters differs by product, you can stabilize the small number of new products with information from existing product lines.

Approach to Analysis and Modeling

The common prior distribution of yield is set at an average of 95%, equivalent to 40 $\mathrm{Beta}(38,2)$, and the posterior distribution by product is calculated.

Check with Python

yield_a0, yield_b0 = 38.0, 2.0
product_df["post_a"] = yield_a0 + product_df["passed"]
product_df["post_b"] = yield_b0 + product_df["produced"] - product_df["passed"]
product_df["hier_yield"] = product_df["post_a"] / (product_df["post_a"] + product_df["post_b"])
product_df["prob_below_94pct"] = stats.beta.cdf(0.94, product_df["post_a"], product_df["post_b"])
product_df

	product	produced	passed	observed_yield	post_a	post_b	hier_yield	prob_below_94pct
0	Product-1	1000	965	0.9650	1003.0	37.0	0.9644	1.4545e-04
1	Product-2	700	664	0.9486	702.0	38.0	0.9486	1.4382e-01
2	Product-3	350	346	0.9886	384.0	6.0	0.9846	3.3491e-06
3	Product-4	160	145	0.9062	183.0	17.0	0.9150	9.0861e-01
4	Product-5	60	58	0.9667	96.0	4.0	0.9600	1.4837e-01
5	Product-6	20	18	0.9000	56.0	4.0	0.9333	5.2445e-01

Reading the results

Zero or a few failures for small quantities of products are not considered permanent product differences, but rather as the level of existing product groups. On the other hand, products with a high probability of not meeting standards are subject to review for process conditions and design tolerances.

No.075: Estimating response rates by customer segment

Meaning in Practice

Customer response rates are replaced by improvement success rates by material supplier segment on the manufacturing floor. Compare the difference in effectiveness by segment by counting the number of trials.

Approach to Analysis and Modeling

A common improvement success rate of 60% is updated from a pre-distribution of 20 cases, and the probability that the success rate exceeds 60% is calculated.

Check with Python

segments = ["Domestic-A", "Domestic-B", "Overseas-A", "Overseas-B"]
trials = np.array([90, 45, 18, 7])
true_response = np.array([0.72, 0.64, 0.55, 0.70])
success = rng.binomial(trials, true_response)
segment_df = pd.DataFrame({"segment": segments, "trials": trials, "success": success})
segment_df["observed"] = segment_df["success"] / segment_df["trials"]
a_seg, b_seg = 12.0, 8.0
segment_df["posterior"] = (a_seg + segment_df["success"]) / (a_seg + b_seg + segment_df["trials"])
segment_df["prob_over_60pct"] = 1 - stats.beta.cdf(0.60, a_seg + segment_df["success"], b_seg + segment_df["trials"] - segment_df["success"])
segment_df

	segment	trials	success	observed	posterior	prob_over_60pct
0	Domestic-A	90	66	0.7333	0.7091	0.9919
1	Domestic-B	45	33	0.7333	0.6923	0.9421
2	Overseas-A	18	7	0.3889	0.5000	0.1080
3	Overseas-B	7	4	0.5714	0.5926	0.4787

Reading the results

Even if a category with fewer cases has a high performance rate, the post-event probability is not necessarily sufficiently high. Company-wide deployment is decided not only by expectations but also by combining minimum accuracy and additional testing costs.

No.076: Stabilizing Estimates for Minority Data Groups

Meaning in Practice

Because there is little data on prototype lines and new products, rankings change significantly every month. By numerically indicating the strength of reduction, you can explain to the site why the correction was made.

Approach to Analysis and Modeling

The weight of performance in the post-average is $w_j=n_j/(n_j+40)$. The remaining $1-w_j$ is the weight of common prior information.

Check with Python

line_df["data_weight"] = line_df["inspected"] / (line_df["inspected"] + alpha0 + beta0)
line_df["shrinkage"] = line_df["observed_rate"] - line_df["hier_rate"]
line_df[["line", "inspected", "observed_rate", "hier_rate", "data_weight", "shrinkage"]]

	line	inspected	observed_rate	hier_rate	data_weight	shrinkage
0	Line-A	1200	0.0200	0.0202	0.9677	-0.0002
1	Line-B	900	0.0211	0.0213	0.9574	-0.0002
2	Line-C	650	0.0338	0.0333	0.9420	0.0005
3	Line-D	400	0.0175	0.0182	0.9091	-0.0007
4	Line-E	180	0.0111	0.0136	0.8182	-0.0025
5	Line-F	80	0.0125	0.0167	0.6667	-0.0042
6	Line-G	35	0.1143	0.0667	0.4667	0.0476
7	Line-H	12	0.0833	0.0385	0.2308	0.0449

fig, ax = plt.subplots(figsize=(8, 4))
ax.scatter(line_df["inspected"], line_df["data_weight"], s=60)
for _, r in line_df.iterrows(): ax.annotate(r["line"], (r["inspected"], r["data_weight"]), xytext=(4, 4), textcoords="offset points")
ax.set_title("Data volume and weight assigned to observed results")
ax.set_xlabel("Number inspected")
ax.set_ylabel("Weight on observed rate")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

Reading the results

The group of 12 tests is mainly estimated from overall trends, while the group of 1,200 is determined almost entirely by its own group data. Judgments for minority groups are not categorized as “normal,” but are linked to additional measurement plans to narrow the estimated range.

No.077: Understanding the Random Slice Model

Meaning in Practice

Even under the same temperature and equipment conditions, there are line-specific standard yields. Random intercepts treat this baseline difference as a deviation from the company-wide distribution.

Approach to Analysis and Modeling

Daily yield y_{ij}=eta_0+u_j+arepsilon_{ij} First, let’s calculate the line average. Normal Methods Using Known Intergroup and Intragroup Dispersions—With a simplified version of the official model, $u_j$ is estimated.

Check with Python

days = 30
intercepts = np.array([0.010, -0.004, 0.006, -0.012, 0.003, 0.015, -0.008, 0.000])
daily = []
for line, u in zip(lines, intercepts):
    for day in range(days):
        daily.append((line, day + 1, 0.95 + u + rng.normal(0, 0.012)))
daily_df = pd.DataFrame(daily, columns=["line", "day", "yield_rate"])
means = daily_df.groupby("line")["yield_rate"].agg(["mean", "count", "std"]).reset_index()
grand = daily_df["yield_rate"].mean()
tau2, sigma2 = 0.010**2, 0.012**2
means["reliability"] = tau2 / (tau2 + sigma2 / means["count"])
means["random_intercept_est"] = means["reliability"] * (means["mean"] - grand)
means

	line	mean	count	std	reliability	random_intercept_est
0	Line-A	0.9613	30	0.0111	0.9542	0.0101
1	Line-B	0.9437	30	0.0121	0.9542	-0.0066
2	Line-C	0.9581	30	0.0079	0.9542	0.0071
3	Line-D	0.9350	30	0.0131	0.9542	-0.0150
4	Line-E	0.9514	30	0.0126	0.9542	0.0006
5	Line-F	0.9620	30	0.0145	0.9542	0.0108
6	Line-G	0.9404	30	0.0106	0.9542	-0.0098
7	Line-H	0.9536	30	0.0149	0.9542	0.0027

Reading the results

If the estimated random intercept is positive, it may still have a higher-than-average reference yield after conditional adjustment. Before concluding that it is a permanent difference on site, we check factors such as equipment generation, team, and material lot that have not been deployed.

No.078: Understanding the Random Coefficient Model

Meaning in Practice

The effects of temperature rise can vary depending on the line. Having only a common coefficient misses local weaknesses, while regression by line picks up noise.

Approach to Analysis and Modeling

y_{ij}=eta_0+u_{0j}+(eta_1+u_{1j})x_{ij}+arepsilon_{ij}, calculate the temperature coefficient for each line, and partially reduce it to the overall coefficient. Here, we provide a simple estimate that clearly explains the mechanism.

Check with Python

slope_true = np.array([-0.0010, -0.0018, -0.0007, -0.0022, -0.0013, -0.0005, -0.0020, -0.0011])
rows = []
for line, b0, b1 in zip(lines, intercepts, slope_true):
    temp = rng.normal(0, 2.5, days)
    y = 0.95 + b0 + b1 * temp + rng.normal(0, 0.008, days)
    rows.extend(zip([line]*days, temp, y))
temp_df = pd.DataFrame(rows, columns=["line", "temp_deviation", "yield_rate"])
est = []
for line, g in temp_df.groupby("line"):
    slope, intercept = np.polyfit(g["temp_deviation"], g["yield_rate"], 1)
    est.append((line, intercept, slope))
slope_df = pd.DataFrame(est, columns=["line", "intercept", "raw_slope"])
global_slope = slope_df["raw_slope"].mean()
slope_df["partial_slope"] = 0.70 * slope_df["raw_slope"] + 0.30 * global_slope
slope_df.sort_values("partial_slope")

	line	intercept	raw_slope	partial_slope
4	Line-E	0.9519	-0.0022	-0.0020
3	Line-D	0.9408	-0.0021	-0.0019
6	Line-G	0.9419	-0.0017	-0.0016
1	Line-B	0.9458	-0.0016	-0.0016
0	Line-A	0.9600	-0.0013	-0.0013
2	Line-C	0.9554	-0.0011	-0.0012
7	Line-H	0.9479	-0.0009	-0.0010
5	Line-F	0.9628	-0.0007	-0.0009

Reading the results

The more negative the coefficient, the more sensitive it is to temperature deviations. Maintenance and air conditioning investments can be prioritized not only by average yield but also by sensitivity to environmental changes.

No.079: Visualizing the Estimation Results of the Hierarchical Bayesian Model

Meaning in Practice

In management and quality meetings, it is necessary to convey not only averages but also the relationship between uncertainty and baseline values in a single sheet.

Approach to Analysis and Modeling

The post-poster average by line and the 95% credit range are used as forest plots, with an administrative attention level of 3%.

Check with Python

viz = line_df.sort_values("hier_rate").reset_index(drop=True)
xerr = np.vstack([viz["hier_rate"] - viz["ci_low"], viz["ci_high"] - viz["hier_rate"]])
fig, ax = plt.subplots(figsize=(9, 5))
ax.errorbar(viz["hier_rate"], viz["line"], xerr=xerr, fmt="o", capsize=4)
ax.axvline(0.03, color="crimson", linestyle="--", label="Attention threshold: 3%")
ax.set_title("Posterior defect rates with 95% credible intervals")
ax.set_xlabel("Defect rate")
ax.set_ylabel("Production line")
ax.grid(axis="x", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

Reading the results

If the entire section exceeds the standard, there is strong justification for improvement intervention; if the section is broad and exceeds the standard, additional testing is reasonable. Instead of ‘highest scores,’ both position and width are read simultaneously.

No.080: Organizing the Criteria for Using Hierarchical Bayes in Practical Data Analysis

Meaning in Practice

Hierarchical Bayes are not always necessary. We decide on adoption by checking the number of groups, data volume per group, commonality, and the cost of misjudgment.

Approach to Analysis and Modeling

Candidate projects are scored based on “hierarchical structure,” “small group,” “inter-group difference,” and “misjudgment cost.” This is not statistical proof, but rather a business check that begins modeling considerations.

Check with Python

use_cases = pd.DataFrame({
    "case": ["Multi-line defect rate", "Single machine stable process", "New product yield", "Supplier improvement test"],
    "hierarchy": [2, 0, 2, 2], "small_groups": [2, 0, 2, 1],
    "between_group_diff": [2, 0, 1, 2], "decision_cost": [2, 1, 2, 2]
})
use_cases["score"] = use_cases.iloc[:, 1:].sum(axis=1)
use_cases["recommendation"] = pd.cut(use_cases["score"], [-1, 2, 5, 8], labels=["Simple model first", "Compare both", "Hierarchical model candidate"])
use_cases.sort_values("score", ascending=False)

	case	hierarchy	small_groups	between_group_diff	decision_cost	score	recommendation
0	Multi-line defect rate	2	2	2	2	8	Hierarchical model candidate
2	New product yield	2	2	1	2	7	Hierarchical model candidate
3	Supplier improvement test	2	1	2	2	7	Hierarchical model candidate
1	Single machine stable process	0	0	0	1	1	Simple model first

Reading the results

Multi-line quality and new product yield are easy candidates to consider, and simple models come first when a single equipment has a sufficient time sequence. At the time of implementation, we compare prediction performance, interpretability, and operational costs with simple models.

Practical Implications Seen Through Target Exercise

1. Judging by probability, not ranking: Use the probability of exceeding standards and credit bands as common indicators for improvement meetings.
2. A minority group is more ‘insufficient’ than ‘good/bad’: Clearly state the value of additional testing.
3. Balancing commonality and individualization: Share company-wide standards while retaining site-specific differences.
4. Look not only at the slices but also at sensitivity differences: Connecting lines vulnerable to environmental changes to preventive conservation.

What is necessary for practical implementation

• Record hierarchical keys such as lines, equipment, varieties, and material lots without missing any omissions.
• Document the basis for prior distributions by dividing them into past periods, expert judgment, and weak information prior distributions.
• Inspecting for biases such as production volume differences, changes in inspection methods, and post-sorting data.
• Conduct sensitivity analysis of pre-distribution data, post-prediction checks, and out-of-period validation.
• Agree on probability thresholds and losses corresponding to “improvement,” “additional measurement,” and “continued monitoring.”
• Include model update frequency, responsible persons, audit logs, and explanatory materials for the field in the operational design.

Conclusion

The value of hierarchical Bayes lies not in the complex model itself but in Fairly handling differences in data volume, and simultaneously incorporating overall knowledge and field differences into decision-making.. It is practical to start with small, explainable implementations such as conjugate models, verify improvements over simple aggregation, and then extend the model.

Consultations for Corporations

At Suri Kobo, we support you with tailored business challenges, from hierarchical design of quality data, Bayesian model PoC, decision-making rule design, to on-site presentation reports.

📩 Contact Us: surikobo.co.jp/contact Please feel free to consult us first.