具有正弦波纹的n阶多项式曲线拟合

Question

我在为某个测量装置的测量误差建模。这就是data的外观：低频多项式上的高频正弦波纹。我的模型也能捕捉到涟漪。在

拟合误差的曲线应为：误差（x）=a0+a1*x+a2*x^2+。。。一个*x^n+Asin（x/lambda）。多项式的阶数是未知的。我的计划是从1-9迭代n并选择一个具有最高的F-value。在

到目前为止，我玩过numpy.polyfit和{}。numpy.polyfit仅适用于多项式，因此虽然我可以生成“最佳拟合”多项式，但无法确定正弦项的参数A和lambda。scipy.optimize.curve_fit如果我已经知道了误差（x）多项式部分的多项式的阶数，那么它的效果会很好。在

有没有一种聪明的方法可以同时使用numpy.polyfit和{}来完成这项工作？或者另一个图书馆功能？

下面是我如何使用numpy.polyfit选择最佳多项式的代码：

def GetErrorPolynomial(X, Y):

    maxFval = 0.0
    for i in range(1, 10):   # i is the order of the polynomial (max order = 9)
        error_func = np.polyfit(X, Y, i)
        error_func = np.poly1d(error_func)

        # F-test (looking for the largest F value)
        numerator = np.sum(np.square(error_func(X) - np.mean(Y))) / i
        denominator = np.sum(np.square(Y - error_func(X))) / (Y.size - i - 1)
        Fval = numerator / denominator

        if Fval > maxFval:
            maxFval = Fval
            maxFvalPolynomial = error_func

    return maxFvalPolynomial

下面是我如何使用curve_fit的代码：

^{pr2}$

它“硬编码”为二次函数，但我想选择“最佳”顺序，就像我上面用np.polyfit所做的那样

Answer 1

使用Scikit-learn Linear Regression

下面是一个代码示例，我使用一个3次多项式执行线性回归，该多项式通过值为1和零导数的点0。您只需将create_vector函数与所需函数相适应。在

from sklearn import linear_model
import numpy as np

def create_vector(x):
    # currently representing a polynom Y = a*X^3 + b*X^2
    x3 = np.power(x, 3)
    x2 = np.power(x, 2)
    X = np.append(x3, x2, axis=1)
    return X 

data_x = [some_data_input]
data_y = [some_data_output]

x = np.array(data_x).reshape(-1, 1)
y_data = np.array(data_y).reshape(-1, 1)-1 # -1 to pass by the point (0,1)

X = create_vector(x)    

regr = linear_model.LinearRegression(fit_intercept=False)
regr.fit(X, y_data)

Answer 2

我终于找到了一种方法来模拟涟漪，并能回答我自己的问题。这个2006 paper对与我的数据集相似的波纹进行曲线拟合。在

首先，我做了一个最小二乘多项式拟合，然后从原始数据中减去这个多项式曲线。这只给我留下了涟漪。应用傅立叶变换，我选择了主频，让我重建正弦波纹。然后我简单地把这些波纹加到我刚开始得到的多项式曲线上。是这样的。在

Answer 3

我从散点图中提取数据进行分析，发现多项式+正弦似乎不是一个最佳模型，因为低阶多项式没有很好地遵循数据的形状，高阶多项式在数据极值处表现出Runge的高曲率现象。（“0.1）可能是（a）和（b）方程中的一个很好的（a）和（0.1）的峰值（b）的搜索结果。在

这是一个图形化的Python曲线拟合器，在文件的顶部我加载了我提取的数据，所以你需要用实际的数据替换它。该算法使用scipy的差分进化遗传算法模块来估计非线性fitter的初始参数值，该算法使用拉丁超立方体算法确保参数空间的彻底搜索，并要求搜索范围。在这里，这些界限是从数据的最大值和最小值获取的。在

从拟合曲线中减去模型预测值后，您将只剩下要建模的正弦分量。我注意到在x=275附近似乎还有一个额外的窄的低振幅峰值。在

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

##########################################################
# load data section
f = open('/home/zunzun/temp/temp.dat')
textData = f.read()
f.close()

xData = []
yData = []
for line in textData.split('\n'):
    if line: # ignore blank lines
        spl = line.split()
        xData.append(float(spl[0]))
        yData.append(float(spl[1]))
xData = numpy.array(xData)
yData = numpy.array(yData)


##########################################################
# model to be fitted
def func(x, a, b, c, offset): # Extreme Valye Peak equation from zunzun.com
    return a * numpy.exp(-1.0 * numpy.exp(-1.0 * ((x-b)/c))-((x-b)/c) + 1.0) + offset


##########################################################
# fitting section

# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)

    minData = min(minX, minY)
    maxData = max(maxX, maxY)

    parameterBounds = []
    parameterBounds.append([minData, maxData]) # search bounds for a
    parameterBounds.append([minData, maxData]) # search bounds for b
    parameterBounds.append([minData, maxData]) # search bounds for c
    parameterBounds.append([minY, maxY]) # search bounds for offset

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)

更新- 如果高频正弦分量是常数（我不知道），那么仅用几个周期对一小部分数据进行建模就足以确定拟合模型正弦波部分的方程和初始参数估计。在这里，我做了这个，结果如下：

根据以下方程式：

^{pr2}$

结合这两个模型使用实际数据，而不是我的散点图提取的数据，应该接近您需要的。在