通过置信区间拟合金属疲劳实验数据的Python代码
时间: 2024-02-17 11:48:33 浏览: 55
基于percols算法的疲劳检测python代码
5星 · 资源好评率100%
以下是通过置信区间拟合金属疲劳实验数据的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()
```
输出结果:
![confidence_interval_fit](https://i.imgur.com/YdRbTJh.png)
阅读全文