Python实战：用NumPy和SciPy玩转最小二乘拟合

传感器数据在采集过程中总会有噪声干扰，实验室测量也难免存在系统误差。当我们需要从这些不完美的数据中提取出有价值的规律时，最小二乘法就成了数据科学家工具箱中的瑞士军刀。本文将带你用Python的NumPy和SciPy库，从零实现最小二乘拟合，并解决实际工程中的校准问题。

1. 最小二乘法的核心思想

最小二乘法的本质是通过最小化误差的平方和来寻找数据的最佳匹配函数。想象一下这样的场景：你有一组实验测量的数据点，想要找到一条直线或曲线，使得这些点到这条曲线的垂直距离之和最小。

为什么选择平方和最小化？

• 平方运算能保证所有误差项为正
• 对大误差的惩罚更严厉（平方放大效应）
• 数学处理方便（可导、有解析解）

在Python中，我们可以用NumPy的线性代数模块高效实现这些计算。先来看一个简单的线性回归例子：

import numpy as np

# 生成带噪声的线性数据
np.random.seed(42)
x = np.linspace(0, 10, 50)
y = 2.5 * x + 1.2 + np.random.normal(scale=1.5, size=len(x))

# 构造设计矩阵
A = np.vstack([x, np.ones(len(x))]).T

# 使用正规方程求解
theta = np.linalg.inv(A.T @ A) @ A.T @ y
slope, intercept = theta

print(f"拟合结果: y = {slope:.2f}x + {intercept:.2f}")

这段代码演示了最小二乘法的矩阵解法核心——正规方程。我们构造设计矩阵A，然后通过矩阵运算直接求出最优参数。

2. NumPy实现正规方程

对于工程师来说，理解正规方程的推导过程比记住最终公式更重要。让我们拆解这个矩阵运算背后的原理。

正规方程的推导过程：

1. 定义误差函数：$S = \sum_{i=1}^n (y_i - (mx_i + b))^2$
2. 用矩阵表示：$S = (Y - X\theta)^T(Y - X\theta)$
3. 对$\theta$求导并令导数为零，得到：$X^TX\theta = X^TY$
4. 解得：$\theta = (X^TX)^{-1}X^TY$

在实际应用中，我们需要注意几个关键点：

• 矩阵条件数：当$X^TX$接近奇异矩阵时，求逆会不稳定
• 数值稳定性：对于病态问题，建议使用QR分解或SVD
• 计算效率：对于大数据集，直接求逆可能很慢

下面是一个更健壮的实现，加入了异常处理：

def linear_least_squares(x, y):
    """
    稳健的线性最小二乘拟合实现
    参数:
        x: 自变量数组
        y: 因变量数组
    返回:
        (斜率, 截距), 拟合优度R²
    """
    A = np.vstack([x, np.ones(len(x))]).T
    
    try:
        # 使用伪逆提高数值稳定性
        theta = np.linalg.pinv(A) @ y
    except np.linalg.LinAlgError:
        print("矩阵奇异，无法求解")
        return None
    
    # 计算R²
    y_pred = A @ theta
    ss_tot = np.sum((y - np.mean(y))**2)
    ss_res = np.sum((y - y_pred)**2)
    r_squared = 1 - (ss_res / ss_tot)
    
    return theta, r_squared

3. SciPy的curve_fit实战

对于更复杂的非线性拟合，SciPy的curve_fit函数提供了便捷的接口。让我们通过一个传感器校准的实际案例来演示它的强大功能。

案例背景： 假设我们有一个温度传感器，其电阻值$R$与温度$T$的关系符合Steinhart-Hart方程： $\frac{1}{T} = A + B\ln(R) + C(\ln(R))^3$

我们采集了以下校准数据：

温度(°C)	电阻(Ω)
0	32000
25	10000
50	4500
75	2500
100	1500

用curve_fit拟合这个模型的代码如下：

from scipy.optimize import curve_fit

def steinhart_hart(r, a, b, c):
    """Steinhart-Hart方程"""
    log_r = np.log(r)
    return 1 / (a + b * log_r + c * log_r**3)

# 实验数据
temps = np.array([0, 25, 50, 75, 100]) + 273.15  # 转换为开尔文
resistances = np.array([32000, 10000, 4500, 2500, 1500])

# 初始参数猜测
initial_guess = [1e-3, 1e-4, 1e-7]

# 执行拟合
params, covariance = curve_fit(steinhart_hart, resistances, 1/temps, p0=initial_guess)

a, b, c = params
print(f"拟合参数: A={a:.3e}, B={b:.3e}, C={c:.3e}")

# 评估拟合质量
predicted_temps = 1 / steinhart_hart(resistances, a, b, c) - 273.15
residuals = (temps - 273.15) - predicted_temps
rmse = np.sqrt(np.mean(residuals**2))
print(f"RMSE: {rmse:.2f}°C")

curve_fit使用技巧：

• 提供合理的初始参数猜测(p0)有助于收敛
• 可以通过bounds参数限制参数范围
• 协方差矩阵的对角线元素给出参数方差估计
• 对于复杂模型，可能需要调整maxfev(最大函数评估次数)

4. 模型评估与残差分析

拟合完成后，我们需要评估模型的质量。常见的评估指标包括：

1. 残差图：检查残差是否随机分布
2. R²值：解释的方差比例
3. RMSE：均方根误差，与原始数据同单位
4. 参数置信区间：从协方差矩阵计算

下面是一个完整的评估流程实现：

import matplotlib.pyplot as plt

def evaluate_fit(x, y, y_pred, params, cov):
    """全面评估拟合结果"""
    residuals = y - y_pred
    
    # 计算关键指标
    ss_res = np.sum(residuals**2)
    ss_tot = np.sum((y - np.mean(y))**2)
    r_squared = 1 - (ss_res / ss_tot)
    rmse = np.sqrt(ss_res / len(y))
    
    # 参数标准差
    perr = np.sqrt(np.diag(cov))
    
    # 绘制结果
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # 拟合曲线与数据点
    ax1.scatter(x, y, label='实际数据')
    ax1.plot(x, y_pred, 'r-', label='拟合曲线')
    ax1.set_title(f"拟合结果 (R²={r_squared:.3f}, RMSE={rmse:.3f})")
    ax1.legend()
    
    # 残差分析
    ax2.scatter(x, residuals)
    ax2.axhline(y=0, color='r', linestyle='--')
    ax2.set_title("残差图")
    
    plt.tight_layout()
    plt.show()
    
    # 打印参数信息
    print("\n参数估计:")
    for i, (param, err) in enumerate(zip(params, perr)):
        print(f"参数 {i}: {param:.5f} ± {err:.5f}")
    
    return r_squared, rmse

解读残差图的关键点：

• 残差应随机分布在零线周围
• 明显的模式可能表明模型选择不当
• 异方差性（残差幅度随x变化）需要关注
• 离群点可能需要特别处理

5. 加权最小二乘法处理异方差数据

当不同数据点的测量精度不同时，我们需要使用加权最小二乘法。常见场景包括：

• 不同仪器采集的数据
• 不同实验条件下的测量
• 随时间变化的测量精度

SciPy的curve_fit支持weights参数来实现加权拟合。下面是一个示例：

# 假设后三个数据点测量更精确
weights = np.array([1, 1, 2, 2, 3])

# 加权拟合
params_weighted, _ = curve_fit(
    steinhart_hart, 
    resistances, 
    1/temps, 
    p0=initial_guess,
    sigma=1/np.sqrt(weights)  # 注意sigma是标准差
)

# 比较结果
print("普通拟合参数:", params)
print("加权拟合参数:", params_weighted)

权重设置的常见策略：

1. 与测量误差的方差成反比
2. 基于数据采集的置信度
3. 对于重复测量，使用测量次数的平方根

6. 处理特殊情况的实用技巧

在实际工程应用中，我们会遇到各种特殊情况。以下是几个常见问题的解决方案：

问题1：矩阵奇异或接近奇异

# 使用SVD分解代替直接求逆
U, s, Vh = np.linalg.svd(A, full_matrices=False)
theta = Vh.T @ np.diag(1/s) @ U.T @ y

问题2：拟合指数衰减模型不收敛

# 对指数模型先取对数转换为线性问题
def exp_model(x, a, b):
    return a * np.exp(b * x)

# 初始参数估计技巧
y_log = np.log(y)
A = np.vstack([x, np.ones(len(x))]).T
b_lin, a_lin = np.linalg.lstsq(A, y_log, rcond=None)[0]
initial_guess = [np.exp(a_lin), b_lin]

问题3：有约束条件的拟合

from scipy.optimize import minimize

def residual(params, x, y):
    a, b = params
    return np.sum((y - (a * x + b))**2)

# 约束b > 0
constraints = {'type': 'ineq', 'fun': lambda params: params[1]}
result = minimize(residual, [1, 0], args=(x, y), constraints=constraints)
a_opt, b_opt = result.x

7. 实际工程案例：多传感器数据融合

在物联网应用中，我们经常需要融合多个传感器的数据。假设有三个温度传感器测量同一环境温度，但精度不同：

# 三个传感器的读数及其标定误差
sensor_data = np.array([
    [25.1, 25.3, 24.8],  # 时间点1
    [26.0, 25.8, 26.5],  # 时间点2
    [24.5, 24.3, 24.7]   # 时间点3
])
sensor_errors = np.array([0.5, 0.3, 0.8])  # 各传感器的标准差

# 计算最优加权平均
weights = 1 / sensor_errors**2
weighted_avg = np.sum(sensor_data * weights, axis=1) / np.sum(weights)

print("融合后的温度估计:", weighted_avg)

这个例子展示了最小二乘思想在传感器融合中的应用。通过给更精确的传感器分配更高权重，我们得到了比简单平均更可靠的结果。

8. 性能优化与大规模数据处理

当处理大规模数据集时，最小二乘计算可能变得昂贵。以下是几种优化策略：

方法1：使用QR分解代替正规方程

Q, R = np.linalg.qr(A)
theta = np.linalg.solve(R, Q.T @ y)

方法2：随机梯度下降(SGD)

def sgd_least_squares(X, y, learning_rate=0.01, epochs=100):
    theta = np.zeros(X.shape[1])
    m = len(y)
    
    for _ in range(epochs):
        for i in range(m):
            xi = X[i:i+1]  # 保持2D形状
            gradient = -2 * xi.T @ (y[i] - xi @ theta)
            theta -= learning_rate * gradient
    
    return theta

方法3：使用稀疏矩阵

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import lsqr

# 当设计矩阵A很大且稀疏时
A_sparse = csr_matrix(A)
theta = lsqr(A_sparse, y)[0]

9. 从二维到多维：多元线性回归

当因变量受多个因素影响时，我们需要扩展模型到多维：

# 假设y取决于x1, x2两个变量
x1 = np.random.rand(100)
x2 = np.random.rand(100)
y = 3 * x1 - 2 * x2 + 1 + np.random.normal(scale=0.1, size=100)

# 构造设计矩阵
A = np.column_stack([x1, x2, np.ones_like(x1)])

# 求解
theta = np.linalg.lstsq(A, y, rcond=None)[0]
print(f"多元回归系数: x1={theta[0]:.2f}, x2={theta[1]:.2f}, 截距={theta[2]:.2f}")

解读结果注意事项：

• 系数大小反映变量重要性，但需考虑变量单位
• 多重共线性会导致系数估计不稳定
• 可以使用标准化处理（z-score）使系数可比

10. 非线性模型选择与比较

选择合适的模型形式是拟合成功的关键。以下是几种常见模型及其适用场景：

1. 多项式模型：

   def polynomial(x, a, b, c):
       return a * x**2 + b * x + c
   

   • 适用场景：光滑曲线，无渐近线

2. 指数衰减模型：

   def exponential_decay(x, a, b, c):
       return a * np.exp(-b * x) + c
   

   • 适用场景：衰减过程，如放射性衰变、电容放电

3. 对数模型：

   def logarithmic(x, a, b):
       return a * np.log(x) + b
   

   • 适用场景：增长递减现象

4. 幂律模型：

   def power_law(x, a, b):
       return a * x**b
   

   • 适用场景：分形、标度现象

模型比较方法：

• AIC/BIC信息准则
• 交叉验证误差
• 残差分析

from sklearn.model_selection import cross_val_score
from sklearn.metrics import make_scorer

def rmse(y_true, y_pred):
    return np.sqrt(np.mean((y_true - y_pred)**2))

# 比较多项式模型次数
results = {}
for degree in [1, 2, 3]:
    model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    scores = cross_val_score(model, x[:, np.newaxis], y, 
                           scoring=make_scorer(rmse), cv=5)
    results[degree] = np.mean(scores)

print("交叉验证RMSE:", results)

11. 稳健回归处理异常值

当数据中存在异常值时，常规最小二乘会受到影响。稳健回归方法通过降低异常值的权重来提高模型鲁棒性。

Huber回归实现：

from sklearn.linear_model import HuberRegressor

# 添加异常值
x_out = np.append(x, [12, 13, 14])
y_out = np.append(y, [50, 60, 70])

# 常规线性回归
lr = LinearRegression().fit(x_out[:, np.newaxis], y_out)

# Huber回归
huber = HuberRegressor(epsilon=1.35).fit(x_out[:, np.newaxis], y_out)

# 比较结果
plt.scatter(x_out, y_out)
plt.plot(x_out, lr.predict(x_out[:, np.newaxis]), 'r-', label='OLS')
plt.plot(x_out, huber.predict(x_out[:, np.newaxis]), 'g--', label='Huber')
plt.legend()

其他稳健回归方法：

• RANSAC（随机抽样一致）
• Theil-Sen估计器
• 最小中值平方回归

12. 贝叶斯视角下的线性回归

贝叶斯方法提供了参数的不确定性估计。使用PyMC3实现：

import pymc3 as pm

with pm.Model() as linear_model:
    # 先验分布
    a = pm.Normal('a', mu=0, sigma=10)
    b = pm.Normal('b', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)
    
    # 线性模型
    mu = a * x + b
    
    # 似然函数
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
    
    # 采样
    trace = pm.sample(2000, tune=1000)

# 查看后验分布
pm.plot_posterior(trace, var_names=['a', 'b'])

贝叶斯方法的优势：

• 自然获得参数的不确定性区间
• 可以融入先验知识
• 适合小样本情况

13. 实时更新参数的递归最小二乘

对于流式数据，我们可以使用递归最小二乘(RLS)来在线更新参数：

class RecursiveLeastSquares:
    def __init__(self, n_features, lambda_=1.0):
        self.n_features = n_features
        self.lambda_ = lambda_  # 遗忘因子
        self.P = np.eye(n_features) / 1e-6  # 逆协方差矩阵
        self.theta = np.zeros(n_features)
    
    def update(self, x, y):
        x = np.asarray(x).reshape(-1, 1)
        k = (self.P @ x) / (self.lambda_ + x.T @ self.P @ x)
        self.theta += k.flatten() * (y - x.T @ self.theta)
        self.P = (self.P - k @ x.T @ self.P) / self.lambda_
        return self.theta

# 使用示例
rls = RecursiveLeastSquares(n_features=2)
for xi, yi in zip(A, y):
    theta = rls.update(xi, yi)

RLS在以下场景特别有用：

• 实时传感器数据处理
• 参数随时间缓慢变化的系统
• 资源受限的嵌入式系统

14. 交叉验证与超参数调优

对于带正则化的回归模型（如岭回归、Lasso），我们需要选择最佳的正则化强度：

from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV

# 定义参数网格
param_grid = {'alpha': np.logspace(-4, 4, 20)}

# 网格搜索
grid = GridSearchCV(Ridge(), param_grid, cv=5, scoring='neg_mean_squared_error')
grid.fit(A, y)

# 最佳参数
print(f"最佳alpha: {grid.best_params_['alpha']:.4f}")
print(f"最佳分数(MSE): {-grid.best_score_:.4f}")

# 绘制正则化路径
alphas = np.logspace(-4, 4, 100)
coefs = []
for a in alphas:
    ridge = Ridge(alpha=a).fit(A, y)
    coefs.append(ridge.coef_)

plt.figure(figsize=(10, 6))
plt.plot(alphas, coefs)
plt.xscale('log')
plt.xlabel('alpha')
plt.ylabel('系数大小')
plt.title('岭回归正则化路径')

选择正则化参数的准则：

• 交叉验证误差最小
• AIC/BIC准则
• 基于领域知识的约束

15. 工业级应用建议

在实际工程中应用最小二乘法时，以下经验值得注意：

1. 数据预处理：

   • 中心化和标准化可以提高数值稳定性
   • 检查并处理缺失值
   • 考虑去除明显异常值

2. 模型验证：

   • 保留部分数据用于最终验证
   • 使用多种评估指标
   • 检查残差是否符合假设

3. 部署考虑：

   • 对于嵌入式系统，考虑定点运算
   • 实现参数更新而非重新拟合
   • 添加合理性检查防止数值溢出

4. 文档记录：

   • 记录使用的模型和假设
   • 保存拟合参数及其不确定性
   • 注明数据来源和采集条件

def save_model_report(model, X, y, filepath):
    """保存模型关键信息到Markdown文件"""
    y_pred = model.predict(X)
    residuals = y - y_pred
    
    with open(filepath, 'w') as f:
        f.write("# 模型报告\n\n")
        f.write(f"**拟合日期**: {datetime.now()}\n\n")
        f.write("## 参数估计\n")
        if hasattr(model, 'coef_'):
            f.write(f"- 斜率: {model.coef_[0]:.4f}\n")
        if hasattr(model, 'intercept_'):
            f.write(f"- 截距: {model.intercept_:.4f}\n")
        
        f.write("\n## 性能指标\n")
        f.write(f"- R²: {r2_score(y, y_pred):.4f}\n")
        f.write(f"- RMSE: {np.sqrt(mean_squared_error(y, y_pred)):.4f}\n")
        
        f.write("\n## 残差分析\n")
        f.write(f"- 均值: {np.mean(residuals):.4f}\n")
        f.write(f"- 标准差: {np.std(residuals):.4f}\n")
        f.write(f"- 最大正残差: {np.max(residuals):.4f}\n")
        f.write(f"- 最大负残差: {np.min(residuals):.4f}\n")

16. 扩展应用：非线性最小二乘

对于更复杂的非线性模型，我们可以使用更高级的优化技术。以拟合正弦波为例：

def sine_model(x, a, f, phi, c):
    """正弦波模型"""
    return a * np.sin(2 * np.pi * f * x + phi) + c

# 生成带噪声的正弦数据
x_sine = np.linspace(0, 2, 100)
y_sine = 2.5 * np.sin(2 * np.pi * 1.3 * x_sine + 0.5) + 1.2
y_sine += np.random.normal(scale=0.3, size=len(x_sine))

# 使用更高级的优化方法
params_sine, _ = curve_fit(
    sine_model, 
    x_sine, 
    y_sine,
    p0=[2, 1, 0, 1],  # 初始猜测很重要
    method='lm',       # Levenberg-Marquardt算法
    maxfev=2000
)

# 可视化结果
plt.scatter(x_sine, y_sine, label='数据')
plt.plot(x_sine, sine_model(x_sine, *params_sine), 'r-', label='拟合')
plt.legend()

非线性拟合的挑战：

• 初始参数猜测很关键
• 可能陷入局部最优
• 计算成本较高
• 参数可能存在相关性

17. 常见问题排查指南

当拟合结果不理想时，可以按照以下步骤排查：

1. 检查数据质量：

   plt.figure(figsize=(12, 3))
   plt.subplot(131)
   plt.hist(y, bins=20)
   plt.title('y分布')
   
   plt.subplot(132)
   plt.plot(x, y, '.')
   plt.title('x-y关系')
   
   plt.subplot(133)
   plt.plot(y_pred, y - y_pred, '.')
   plt.title('残差图')
   

2. 验证模型假设：

   • 线性性：预测值与残差应无趋势
   • 同方差性：残差方差应恒定
   • 正态性：残差应近似正态分布

3. 数值问题排查：

   • 检查矩阵条件数：np.linalg.cond(A.T @ A)
   • 查看参数尺度：print("参数尺度:", params)
   • 尝试不同的求解方法

4. 模型复杂度分析：

   • 学习曲线：观察训练/验证误差随数据量的变化
   • 特征重要性：分析各变量的贡献度

18. 与机器学习管道的集成

最小二乘回归可以无缝集成到现代机器学习工作流中：

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

# 创建管道
pipe = Pipeline([
    ('poly', PolynomialFeatures(degree=2)),  # 多项式特征
    ('scaler', StandardScaler()),            # 标准化
    ('reg', LinearRegression())              # 线性回归
])

# 数据分割
X_train, X_test, y_train, y_test = train_test_split(
    x[:, np.newaxis], y, test_size=0.2, random_state=42)

# 训练与评估
pipe.fit(X_train, y_train)
score = pipe.score(X_test, y_test)
print(f"测试集R²: {score:.3f}")

# 保存管道
from joblib import dump
dump(pipe, 'regression_pipeline.joblib')

集成时的注意事项：

• 确保预处理只在训练数据上拟合
• 监控数据分布变化
• 考虑使用更复杂的基函数（如样条）

19. 可视化技巧与结果展示

优秀的可视化能更有效地传达拟合结果：

def plot_regression_results(x, y, y_pred, ax=None):
    if ax is None:
        ax = plt.gca()
    
    # 绘制数据和拟合线
    ax.scatter(x, y, color='blue', alpha=0.5, label='数据')
    ax.plot(x, y_pred, color='red', linewidth=2, label='拟合')
    
    # 添加置信区间
    residuals = y - y_pred
    sigma = np.std(residuals)
    ax.fill_between(x, 
                   y_pred - 2*sigma, 
                   y_pred + 2*sigma, 
                   color='pink', alpha=0.3, label='95%区间')
    
    # 添加统计信息
    r2 = r2_score(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    textstr = f'$R^2 = {r2:.3f}$\n$RMSE = {rmse:.3f}$'
    props = dict(boxstyle='round', facecolor='white', alpha=0.5)
    ax.text(0.05, 0.95, textstr, transform=ax.transAxes,
            verticalalignment='top', bbox=props)
    
    ax.legend()
    return ax

# 使用示例
fig, ax = plt.subplots(figsize=(8, 6))
plot_regression_results(x, y, A @ theta, ax=ax)
ax.set_title('回归结果可视化')
ax.set_xlabel('自变量')
ax.set_ylabel('因变量')

高级可视化技巧：

• 部分依赖图展示变量影响
• 残差QQ图检验正态性
• 学习曲线诊断偏差-方差权衡

20. 从理论到实践：完整项目示例

让我们通过一个完整的温度传感器校准项目，串联所有知识点：

# 1. 数据采集
raw_data = pd.read_csv('sensor_calibration.csv')
temps = raw_data['reference_temp'].values
readings = raw_data['sensor_reading'].values

# 2. 探索性分析
plt.scatter(readings, temps)
plt.xlabel('传感器读数 (Ω)')
plt.ylabel('参考温度 (°C)')

# 3. 模型选择 - 使用Steinhart-Hart方程
def steinhart_hart(r, a, b, c):
    log_r = np.log(r)
    return 1 / (a + b * log_r + c * log_r**3)

# 4. 拟合模型
p0 = [1e-3, 1e-4, 1e-7]  # 基于物理知识的初始猜测
params, pcov = curve_fit(steinhart_hart, readings, temps + 273.15, p0=p0)

# 5. 评估拟合
temps_kelvin = temps + 273.15
predicted = steinhart_hart(readings, *params)
residuals = temps_kelvin - predicted

print(f"平均绝对误差: {np.mean(np.abs(residuals)):.2f}K")
print(f"最大误差: {np.max(np.abs(residuals)):.2f}K")

# 6. 保存校准参数
calibration = {
    'date': datetime.now().isoformat(),
    'model': 'Steinhart-Hart',
    'parameters': list(params),
    'covariance': pcov.tolist(),
    'sensor_id': 'THS-001'
}

with open('calibration.json', 'w') as f:
    json.dump(calibration, f, indent=2)

# 7. 实现校准函数
def calibrate_temperature(reading, params):
    """将传感器读数转换为温度(°C)"""
    temp_k = steinhart_hart(reading, *params)
    return temp_k - 273.15

# 8. 验证新数据
new_readings = np.array([8500, 12000, 3000])
calibrated_temps = calibrate_temperature(new_readings, params)
print("校准后的温度:", calibrated_temps)

这个完整流程展示了从数据采集到模型部署的全过程，涵盖了最小二乘法的实际工程应用。