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

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

TL;博士

您的数据集距离线性太远，无法通过稳健回归进行独特处理。在能够回归模数之前，您肯定需要数据集预处理。

稳健回归

首先我们处理您的数据：

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'>}

强异常值不影响回归的情况：

但是如果我们放大，我们至少会检测到两种不同的行为：

结论

• 稳健回归将帮助您过滤掉一些强异常值，而无需手动处理它们
• 您的数据集在某种程度上是非线性的：
  • 第一个数据集表现出负曲率
  • 第二个数据集至少包含两个在分析之前需要拆分的行为，或者如果您希望将其自动化，请在拟合之前添加聚类
• 您需要解决是否必须拟合截距，通常应变-应力曲线经过原点（无应变，无应力）。