R/d2及S/C4估计总体标准差，比较其CPK及规格限概率的差异

news2024/11/28 15:43:06

R/d2 和 S/C4 是用于估计总体标准差的无偏估计方法，通常用于控制图中。这些估计方法的主要目的是通过样本数据来估计总体标准差，以便监测过程的稳定性和变异性，而不需要收集整个总体的数据。

具体来说：

R图中的 R/d2 和 S图中的 S/C4 都是无偏估计，其中 d2 和 C4 是常数，用于修正范围和标准差以获得更接近总体标准差的估计值。

import numpy as np
import scipy.stats as stats
from scipy.stats import norm

# 输入数据
data = [
    [6.4, 7.0, 6.4, 6.4, 7.1],
    [6.8, 6.4, 6.4, 6.3, 6.5],
    [6.3, 7.1, 6.5, 6.4, 7.0],
    [6.1, 6.8, 5.9, 5.8, 6.0],
    [6.4, 6.9, 6.8, 6.5, 6.9],
    [6.6, 6.0, 6.1, 6.2, 5.9],
    [6.3, 6.9, 6.6, 6.2, 6.8],
    [6.4, 5.6, 6.2, 6.0, 5.8],
    [6.3, 6.7, 6.6, 6.4, 6.3],
    [6.7, 5.9, 5.8, 6.3, 6.2],
    [6.6, 7.0, 6.5, 6.4, 7.1],
    [6.8, 6.2, 6.5, 6.2, 5.8]
]


# 3. Shapiro-Wilk检验
shapiro_stat, shapiro_p = stats.shapiro(data)
print("\nShapiro-Wilk检验统计值:", shapiro_stat)
print("Shapiro-Wilk检验p-value:", shapiro_p)
if shapiro_p > 0.05:
    print("数据可能来自正态分布")
else:
    print("数据可能不来自正态分布")


# 计算整体标准差
population_std = np.std(data)

# 计算样本标准差
sample_std = np.std(data, ddof=1)  # 使用ddof参数来指定自由度

print("\ndata整体标准差:", population_std)
print("data样本标准差:", sample_std)

# 计算R/d2估计的总体标准差
r_values = [max(subgroup) - min(subgroup) for subgroup in data]
d2 = 2.326  # 从表格或标准文献中查找
r_bar = np.mean(r_values)
sigma_r = r_bar / d2

# 计算S/C4估计的总体标准差
s_values = [np.std(subgroup, ddof=1) for subgroup in data]
C4 = 0.94  # 从表格或标准文献中查找
s_bar = np.mean(s_values)
sigma_s = s_bar / C4

# 输入过程上下限
upper_spec_limit = 7
lower_spec_limit = 5.5

# 计算CPK
cpk_r = min((upper_spec_limit - np.mean(data)) / (3 * sigma_r), (np.mean(data) - lower_spec_limit) / (3 * sigma_r))
cpk_s = min((upper_spec_limit - np.mean(data)) / (3 * sigma_s), (np.mean(data) - lower_spec_limit) / (3 * sigma_s))

print("\n通过R/d2估计的总体标准差 (σ):", sigma_r)
print("通过S/C4估计的总体标准差 (σ):", sigma_s)
print("R/d2法计算的CPK:", cpk_r)
print("S/C4法计算的CPK:", cpk_s)



# 计算标准分数
z_upper_r = (upper_spec_limit - np.mean(data)) / sigma_r  # 使用R/d2估计的σ
z_lower_r = (lower_spec_limit - np.mean(data)) / sigma_r  # 使用R/d2估计的σ

z_upper_s = (upper_spec_limit - np.mean(data)) / sigma_s  # 使用S/C4估计的σ
z_lower_s = (lower_spec_limit - np.mean(data)) / sigma_s  # 使用S/C4估计的σ

# 计算上限以上的概率（使用R/d2估计的σ）
probability_above_upper_r = 1 - norm.cdf(z_upper_r)

# 计算下限以下的概率（使用R/d2估计的σ）
probability_below_lower_r = norm.cdf(z_lower_r)

# 计算在规格限内的概率（使用R/d2估计的σ）
probability_within_spec_r = 1 - probability_above_upper_r - probability_below_lower_r

# 计算上限以上的概率（使用S/C4估计的σ）
probability_above_upper_s = 1 - norm.cdf(z_upper_s)

# 计算下限以下的概率（使用S/C4估计的σ）
probability_below_lower_s = norm.cdf(z_lower_s)

# 计算在规格限内的概率（使用S/C4估计的σ）
probability_within_spec_s = 1 - probability_above_upper_s - probability_below_lower_s

print("\n使用R/d2法估计的概率(规格上限以上):", probability_above_upper_r)
print("使用R/d2法估计的概率(规格下限以下):", probability_below_lower_r)
print("使用R/d2法估计的概率(在规格限内):", probability_within_spec_r)
print("\n使用S/C4法估计的概率(规格上限以上):", probability_above_upper_s)
print("使用S/C4法估计的概率(规格下限以下):", probability_below_lower_s)
print("使用S/C4法估计的概率(在规格限内):", probability_within_spec_s)

Shapiro-Wilk检验统计值: 0.9730015993118286
Shapiro-Wilk检验p-value: 0.20416395366191864
数据可能来自正态分布

data整体标准差: 0.3711094477673968
data样本标准差: 0.37424122858811426

通过R/d2估计的总体标准差 (σ): 0.3116938950988822
通过S/C4估计的总体标准差 (σ): 0.3245199375603761
R/d2法计算的CPK: 0.6238314176245208
S/C4法计算的CPK: 0.5991756497496196

使用R/d2法估计的概率(规格上限以上): 0.030638302860941935
使用R/d2法估计的概率(规格下限以下): 0.0016361976108266597
使用R/d2法估计的概率(在规格限内): 0.9677254995282314

使用S/C4法估计的概率(规格上限以上): 0.03612600187890791
使用S/C4法估计的概率(规格下限以下): 0.0023663543889196424
使用S/C4法估计的概率(在规格限内): 0.9615076437321725
[Finished in 4.9s]

-------------------------
绘制X-Bar和R图，及数据集直方图概率密度曲线

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 输入数据
data = [
    [6.4, 7.0, 6.4, 6.4, 7.1],
    [6.8, 6.4, 6.4, 6.3, 6.5],
    [6.3, 7.1, 6.5, 6.4, 7.0],
    [6.1, 6.8, 5.9, 5.8, 6.0],
    [6.4, 6.9, 6.8, 6.5, 6.9],
    [6.6, 6.0, 6.1, 6.2, 5.9],
    [6.3, 6.9, 6.6, 6.2, 6.8],
    [6.4, 5.6, 6.2, 6.0, 5.8],
    [6.3, 6.7, 6.6, 6.4, 6.3],
    [6.7, 5.9, 5.8, 6.3, 6.2],
    [6.6, 7.0, 6.5, 6.4, 7.1],
    [6.8, 6.2, 6.5, 6.2, 5.8]
]

# 输入过程上下限
upper_spec_limit = 7
lower_spec_limit = 5.5

# 控制图参数
A2 = 0.577
D4 = 2.113
D3 = 0
d2 = 2.326  # 从表格或标准文献中查找



# 计算X-Bar和R
x_bar = np.mean(data, axis=1)
r_values = np.ptp(data, axis=1)

# 计算X-Bar和R的平均值
x_bar_bar = np.mean(x_bar)
r_bar = np.mean(r_values)
# r_values = [max(subgroup) - min(subgroup) for subgroup in data] #极差均值
# x_double_bar = np.mean([np.mean(subgroup) for subgroup in data]) #x_bar_bar中心线均值


# 将数据展开为一维数组,用于画data数据集直方图
data_flat = [item for sublist in data for item in sublist]
# 计算整体标准差
# population_std = np.std(data_flat)
population_std = r_bar / d2


# 计算UCL和LCL (X-Bar)
UCL_x_bar = x_bar_bar + A2 * r_bar
LCL_x_bar = x_bar_bar - A2 * r_bar

# 计算UCL和LCL (R)
UCL_r = D4 * r_bar
LCL_r = D3 * r_bar

# 绘制X-Bar控制图
plt.figure(figsize=(6, 6))
plt.subplot(3, 1, 1)

plt.plot(x_bar, 'o-', label='X-Bar')
plt.axhline(x_bar_bar, color='r', linestyle='--', label='X-Bar̄')
plt.axhline(UCL_x_bar, color='g', linestyle='--', label='UCL(X-Bar)')
plt.axhline(LCL_x_bar, color='g', linestyle='--', label='LCL(X-Bar)')
plt.ylabel('X-Bar')
# plt.legend()

plt.subplot(3, 1, 2)
plt.plot(r_values, 'o-', color='b', label='R')
plt.axhline(r_bar, color='r', linestyle='--', label='R̄')
plt.axhline(UCL_r, color='g', linestyle='--', label='UCL(R)')
plt.axhline(LCL_r, color='g', linestyle='--', label='LCL(R)')
# plt.xlabel('Sample')
plt.ylabel('R')
# plt.legend()
# plt.title('X-Bar-R')

# 绘制整体数据集的直方图并叠加概率密度曲线
plt.subplot(3, 1, 3)
plt.hist(data_flat, bins=12, density=True, alpha=0.6, color='b', label='Histogram')
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, np.mean(data_flat), population_std)
plt.plot(x, p, 'k', linewidth=2, label='PDF (Population)')
plt.axvline(x_bar_bar, color='r', linestyle='--', label='X-Bar̄')
plt.axvline(UCL_x_bar, color='g', linestyle='--', label='UCL(X-Bar)')
plt.axvline(LCL_x_bar, color='g', linestyle='--', label='LCL(X-Bar)')
plt.axvline(upper_spec_limit, color='b', linestyle='-', label='USL')
plt.axvline(lower_spec_limit, color='b', linestyle='-', label='LSL')
plt.xlabel('Value')
plt.ylabel('Probability')
# plt.legend()
plt.title('Histogram (Population)')
plt.tight_layout()
plt.show()


print("np.std(data_flat)估计总体标准差",np.std(data_flat))
print("r_bar/d2估计总体标准差",r_bar / d2)