通过置信区间拟合金属疲劳实验数据的Python代码

以下是通过置信区间拟合金属疲劳实验数据的Python代码：

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

# 金属疲劳实验数据
x = np.array([8, 12, 16, 20, 24, 28])
y = np.array([13.4, 10.8, 8.4, 6.4, 5.2, 4.2])

# 求解拟合直线的斜率和截距
n = len(x)
x_mean = np.mean(x)
y_mean = np.mean(y)
xy_mean = np.mean(x*y)
x2_mean = np.mean(x**2)
b = (xy_mean - x_mean*y_mean) / (x2_mean - x_mean**2)
a = y_mean - b*x_mean

# 计算拟合直线的预测值和标准误差
y_pred = a + b*x
s = np.sqrt(np.sum((y - y_pred)**2) / (n - 2))
se_b = s / np.sqrt(np.sum((x - x_mean)**2))
se_a = s * np.sqrt((1/n) + (x_mean**2) / np.sum((x - x_mean)**2))

# 计算置信区间
alpha = 0.05
t_val = t.ppf(1 - alpha/2, n-2)
ci_b = [b - t_val*se_b, b + t_val*se_b]
ci_a = [a - t_val*se_a, a + t_val*se_a]

# 绘制数据点和拟合直线
plt.scatter(x, y)
plt.plot(x, y_pred, color='r')

# 绘制置信区间
x_min = np.min(x)
x_max = np.max(x)
y_min = ci_a[0] + ci_b[0]*x_min
y_max = ci_a[1] + ci_b[1]*x_max
plt.plot([x_min, x_max], [y_min, y_max], linestyle='--', color='g')
plt.fill_between([x_min, x_max], [y_min, y_max], alpha=0.1, color='g')

plt.show()

输出结果：