线性拟合前异常数据剔除与平滑的最佳实践及Python代码

我需要通过应力应变曲线的线性拟合计算模量，但采集的压力数据离散度较高，导致拟合效果不理想。计划通过两种方式优化：

• 剔除部分极大/极小值
• 拟合前对压力数据做平滑处理（如滑动平均）

另外，第二组数据的压力文件第180行存在异常点，需要从拟合中排除。

现有拟合代码 #

import sys
import numpy as np
from scipy.optimize import curve_fit
from pathlib import Path
import os


def f(x, a, b):
    return a + b * x

def get_modulus(BoxValue, PressureValue, strain_cutoff = 0.012):
    strain_all = np.array(BoxValue / BoxValue[0] - 1)
    stress_all = np.array(- PressureValue / 10000)
    strain = strain_all[strain_all <= strain_cutoff]
    stress = stress_all[strain_all <= strain_cutoff]
    popt, pcov = curve_fit(f, strain, stress)
    return popt[1], np.sqrt(np.diag(pcov))[1]

# 加载数据示例
BoxValue = np.loadtxt("box-xx_0water_150peg.dat")[:,1]
PressureValue = np.loadtxt(f"pres-xx_0water_150peg")[:,1]
modulus, modulus_std = get_modulus(BoxValue, PressureValue)

数据处理最佳实践与改进代码 #

1. 异常点处理 #

手动指定剔除已知异常点

针对第二组数据第180行的异常点（Python索引为179），加载数据时直接过滤，注意同步剔除对应box数据以保证序列对应：

# 加载第二组压力数据并剔除第180行（索引179）
PressureValue = np.loadtxt("第二组压力文件名")[:,1]
PressureValue = np.delete(PressureValue, 179)

# 同步剔除box数据对应行
BoxValue = np.loadtxt("第二组box文件名")[:,1]
BoxValue = np.delete(BoxValue, 179)

自动剔除极值（通用方案）

用**四分位距(IQR)**法自动识别并剔除离群值，适合未知异常点的场景：

def remove_outliers(data, threshold=1.5):
    q1 = np.percentile(data, 25)
    q3 = np.percentile(data, 75)
    iqr = q3 - q1
    lower_bound = q1 - threshold * iqr
    upper_bound = q3 + threshold * iqr
    mask = (data >= lower_bound) & (data <= upper_bound)
    return mask

# 在拟合函数中加入自动剔除逻辑
def get_modulus(BoxValue, PressureValue, strain_cutoff = 0.012, remove_outlier=True):
    strain_all = np.array(BoxValue / BoxValue[0] - 1)
    stress_all = np.array(- PressureValue / 10000)
    
    # 先筛选应变阈值内的数据
    strain_mask = strain_all <= strain_cutoff
    strain = strain_all[strain_mask]
    stress = stress_all[strain_mask]
    
    # 剔除异常值
    if remove_outlier:
        outlier_mask = remove_outliers(stress)
        strain = strain[outlier_mask]
        stress = stress[outlier_mask]
    
    popt, pcov = curve_fit(f, strain, stress)
    return popt[1], np.sqrt(np.diag(pcov))[1]

2. 数据平滑处理（滑动平均） #

用滑动平均降低数据噪声，平滑后再拟合，窗口大小可根据数据离散程度调整：

def moving_average(data, window_size=5):
    # 简单滑动平均，保持输出长度与输入一致
    weights = np.ones(window_size) / window_size
    smoothed = np.convolve(data, weights, mode='same')
    return smoothed

# 修改拟合函数加入平滑逻辑
def get_modulus(BoxValue, PressureValue, strain_cutoff = 0.012, smooth_window=5):
    strain_all = np.array(BoxValue / BoxValue[0] - 1)
    stress_all = np.array(- PressureValue / 10000)
    
    strain = strain_all[strain_all <= strain_cutoff]
    stress = stress_all[strain_all <= strain_cutoff]
    
    # 平滑应力数据
    stress_smoothed = moving_average(stress, smooth_window)
    
    popt, pcov = curve_fit(f, strain, stress_smoothed)
    return popt[1], np.sqrt(np.diag(pcov))[1]

3. 整合方案：先剔除异常值再平滑 #

结合两种方法进一步提升拟合稳定性，同时增加R²指标评估拟合效果：

def get_modulus_optimized(BoxValue, PressureValue, strain_cutoff = 0.012, 
                          outlier_threshold=1.5, smooth_window=5):
    strain_all = np.array(BoxValue / BoxValue[0] - 1)
    stress_all = np.array(- PressureValue / 10000)
    
    # 筛选应变阈值内的数据
    strain_mask = strain_all <= strain_cutoff
    strain = strain_all[strain_mask]
    stress = stress_all[strain_mask]
    
    # 剔除异常值
    outlier_mask = remove_outliers(stress, outlier_threshold)
    strain = strain[outlier_mask]
    stress = stress[outlier_mask]
    
    # 平滑处理
    stress_smoothed = moving_average(stress, smooth_window)
    
    popt, pcov = curve_fit(f, strain, stress_smoothed)
    # 计算R²评估拟合优度
    r_squared = 1 - (np.sum((stress_smoothed - f(strain, *popt))**2) / 
                     np.sum((stress_smoothed - np.mean(stress_smoothed))**2))
    return popt[1], np.sqrt(np.diag(pcov))[1], r_squared

额外建议 #

• 可视化先行：用matplotlib绘制原始应力应变曲线，直观观察异常点和离散程度，帮助确定阈值和窗口大小：

  import matplotlib.pyplot as plt
  plt.scatter(strain, stress, s=5, label='Raw Data')
  plt.plot(strain, f(strain, *popt), 'r-', label='Fitted Line')
  plt.xlabel('Strain')
  plt.ylabel('Stress')
  plt.legend()
  plt.show()
  

• 加权拟合：如果数据离散程度不均匀，可通过curve_fit的sigma参数给离散小的数据更高权重，提升拟合精度。
• 验证稳定性：尝试不同的异常值阈值、平滑窗口大小，对比模量结果的一致性，选择最优参数。

内容的提问来源于stack exchange，提问作者FreeAir