import numpy as np
import matplotlib.pyplot as plt

# ==================== 数据录入 ====================
# 20个样本，每个样本5个观测值（单位：mm）
data = [
    [29.970, 30.017, 29.898, 29.937, 29.992],
    [29.947, 30.013, 29.993, 29.997, 30.079],
    [30.050, 30.031, 29.999, 29.963, 30.045],
    [30.064, 30.061, 30.016, 30.041, 30.006],
    [29.948, 30.009, 29.962, 29.990, 29.979],
    [30.016, 29.989, 29.939, 29.981, 30.017],
    [29.946, 30.057, 29.992, 29.973, 29.955],
    [29.981, 30.023, 29.992, 29.992, 29.941],
    [30.043, 29.985, 30.014, 29.986, 30.000],
    [30.013, 30.046, 30.096, 29.975, 30.019],
    [30.043, 30.003, 30.062, 30.025, 30.023],
    [29.994, 30.056, 30.033, 30.011, 29.948],
    [29.995, 30.014, 30.018, 29.966, 30.000],
    [30.018, 29.982, 30.028, 30.029, 30.044],
    [30.018, 29.994, 29.995, 30.029, 30.034],
    [30.025, 29.951, 30.038, 30.009, 30.003],
    [30.048, 30.046, 29.995, 30.053, 30.043],
    [30.030, 30.054, 29.997, 29.993, 30.010],
    [29.991, 30.001, 30.041, 30.036, 29.992],
    [30.022, 30.021, 30.022, 30.008, 30.019]
]

# ==================== 计算统计量 ====================
n = 5                    # 子组大小
k = len(data)            # 子组数
means = [np.mean(sample) for sample in data]
ranges = [np.max(sample) - np.min(sample) for sample in data]

x_double_bar = np.mean(means)        # 总均值
r_bar = np.mean(ranges)              # 平均极差

# 控制图系数（n=5）
A2 = 0.577
D3 = 0
D4 = 2.114

# 均值图控制限
UCL_x = x_double_bar + A2 * r_bar
LCL_x = x_double_bar - A2 * r_bar

# 极差图控制限
UCL_r = D4 * r_bar
LCL_r = D3 * r_bar

# ==================== 判断失控点 ====================
out_of_control_x = [i+1 for i, m in enumerate(means) if m > UCL_x or m < LCL_x]
out_of_control_r = [i+1 for i, r in enumerate(ranges) if r > UCL_r]

print("====== 控制图参数 ======")
print(f"总均值 (X̄̄): {x_double_bar:.4f} mm")
print(f"平均极差 (R̄): {r_bar:.4f} mm")
print(f"均值图控制限: UCL={UCL_x:.4f}, LCL={LCL_x:.4f}")
print(f"极差图控制限: UCL={UCL_r:.4f}, LCL={LCL_r:.4f}\n")

if out_of_control_x:
    print(f"⚠️ 均值图失控点（样本编号）: {out_of_control_x}")
else:
    print("✅ 均值图所有点均在控制限内，过程均值受控。")

if out_of_control_r:
    print(f"⚠️ 极差图失控点（样本编号）: {out_of_control_r}")
else:
    print("✅ 极差图所有点均在控制限内，过程离散程度受控。")

# ==================== 绘图 ====================
plt.rcParams['font.sans-serif'] = ['SimHei']   # 用于正常显示中文
plt.rcParams['axes.unicode_minus'] = False     # 正常显示负号

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# ---- 均值图 ----
ax1.plot(range(1, k+1), means, 'bo-', label='样本均值')
ax1.axhline(y=x_double_bar, color='green', linestyle='-', label='CL = {:.4f}'.format(x_double_bar))
ax1.axhline(y=UCL_x, color='red', linestyle='--', label='UCL = {:.4f}'.format(UCL_x))
ax1.axhline(y=LCL_x, color='red', linestyle='--', label='LCL = {:.4f}'.format(LCL_x))
# 标记失控点
for i in out_of_control_x:
    ax1.plot(i, means[i-1], 'rs', markersize=8, label='失控点' if i==out_of_control_x[0] else "")
ax1.set_title('X̄ 控制图 (均值图)')
ax1.set_xlabel('样本组号')
ax1.set_ylabel('均值 (mm)')
ax1.legend()
ax1.grid(True)

# ---- 极差图 ----
ax2.plot(range(1, k+1), ranges, 'bo-', label='样本极差')
ax2.axhline(y=r_bar, color='green', linestyle='-', label='CL = {:.4f}'.format(r_bar))
ax2.axhline(y=UCL_r, color='red', linestyle='--', label='UCL = {:.4f}'.format(UCL_r))
if LCL_r > 0:
    ax2.axhline(y=LCL_r, color='red', linestyle='--', label='LCL = {:.4f}'.format(LCL_r))
for i in out_of_control_r:
    ax2.plot(i, ranges[i-1], 'rs', markersize=8, label='失控点' if i==out_of_control_r[0] else "")
ax2.set_title('R 控制图 (极差图)')
ax2.set_xlabel('样本组号')
ax2.set_ylabel('极差 (mm)')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.show()

# ==================== 工序能力分析（附加） ====================
# 规格限
USL = 30.125
LSL = 29.875
# 估计过程标准差 (d2 for n=5 is 2.326)
d2 = 2.326
sigma_hat = r_bar / d2
Cp = (USL - LSL) / (6 * sigma_hat)
Cpk = min((USL - x_double_bar) / (3 * sigma_hat), (x_double_bar - LSL) / (3 * sigma_hat))

print("\n====== 工序能力分析 ======")
print(f"规格上限 USL = {USL} mm, 规格下限 LSL = {LSL} mm")
print(f"估计标准差 σ̂ = {sigma_hat:.4f} mm")
print(f"Cp = {Cp:.4f}")
print(f"Cpk = {Cpk:.4f}")
if Cpk < 1.33:
    print("⚠️ 工序能力不足 (Cpk < 1.33)，需要改进过程或调整均值。")
else:
    print("✅ 工序能力充足。")