如何排除异常数据点并平滑数据以进行线性拟合

Question

我想通过对应变-应力曲线进行线性拟合来计算模量。不过，由于获得的压力数据波动很大，有时候拟合结果并不好。

我认为有两个方法可以改善这种情况。首先，需要忽略一些极大或极小的值；其次，可能在拟合之前对压力数据进行一些处理会有帮助。例如，使用平滑处理或移动平均可能会有所帮助。

下面是我用来进行线性拟合的Python代码。

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)
   # modulus = popt[1]
   # modulus_std_dev = np.sqrt(np.diag(pcov))[1]
   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)

这里可以下载两组数据。

第一组：盒子数据，和 压力数据

第二组：盒子数据，和 压力数据

第二组的压力数据在第180行有一个异常点，我认为应该在拟合时将其排除。

你能告诉我最佳的处理方法吗？并提供相应的代码吗？

欢迎任何建议或评论。

Answer 1

总结

你的数据集和线性关系相差太远，不能仅靠稳健回归来处理。你需要先对数据集进行预处理，才能进行回归分析。

稳健回归

首先，我们需要处理你的数据：

def prepare(suffix):
    # Load:
    x = pd.read_csv("box" + suffix, sep="\t", header=None, names=["id", "x"])
    y = pd.read_csv("pres" + suffix, sep="\t", header=None, names=["id", "y"])
    # Merge:
    data = x.merge(y, on=["id"])
    # Post process:
    data["strain"] = data["x"] / data["x"][0] - 1.
    data["stress"] = - data["y"] / 10_000.
    return data

然后我们对处理后的数据进行稳健线性回归：

def analyse(data):
    # Regress:
    X, y = data["strain"].values.reshape(-1, 1), data["stress"].values
    regressor = TheilSenRegressor()
    regressor.fit(X, y)
    # Predict:
    yhat = regressor.predict(X)
    score = r2_score(y, yhat)
    # Render:
    fig, axe = plt.subplots()
    axe.scatter(X, y)
    axe.plot(X, yhat, color="orange")
    axe.set_title("Regression")
    axe.set_xlabel("Strain")
    axe.set_ylabel("Stress")
    axe.grid()
    return {
        "slope": regressor.coef_,
        "intercept": regressor.intercept_,
        "score": score,
        "axe": axe
    }

第一个数据集的结果（可能根本不是线性的）：

df1 = prepare("-xx_0water_150peg.dat")
sol1 = analyse(df1)
#{'slope': array([2.22561879]),
# 'intercept': 0.025513888992521497,
# 'score': 0.7902902090920214,
# 'axe': <AxesSubplot:title={'center':'Regression'}, xlabel='Strain', ylabel='Stress'>}

第二个数据集的结果（多种情况）：

df2 = prepare("-yy_200water_150peg.dat")
sol2 = analyse(df2)
#{'slope': array([29.5666213]),
# 'intercept': 0.0008528267952350494,
# 'score': 0.007683701094953643,
# 'axe': <AxesSubplot:title={'center':'Regression'}, xlabel='Strain', ylabel='Stress'>}

在这里，强烈的异常值不会影响回归结果：

但是如果我们放大查看，就能发现至少有两种不同的行为：

结论

• 稳健回归可以帮助你过滤掉一些强烈的异常值，而不需要手动处理它们。
• 你的数据集在某种程度上是非线性的：
  • 第一个数据集表现出负曲率。
  • 第二个数据集包含至少两种行为，在分析之前需要进行拆分，或者如果你想自动化处理，可以在拟合之前添加聚类。
• 你需要考虑是否需要拟合截距，通常应变-应力曲线会经过原点（没有应变就没有应力）。