Ajustement de courbe avec polynôme d'ordre n ayant des ondulations sinusoïdales

Utilisez Régression linéaire Scikit-learn Voici un exemple de code que j'ai utilisé pour effectuer une régression linéaire avec un polynôme de degré 3 qui passe par le point 0 avec la valeur 1 et la dérivée nulle. Il vous suffit d'adapter la fonction create_vector avec la fonction que vous souhaitez.

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)

J'ai extrait des données du nuage de points pour analyse et j'ai trouvé qu'un polynôme + sinus ne semblait pas être un modèle optimal, car les polynômes d'ordre inférieur ne suivaient pas très bien la forme des données et les polynômes d'ordre supérieur présentaient le phénomène de courbure élevée aux données extrêmes. J'ai effectué une recherche d'équation pour trouver à quoi l'onde sinusoïdale haute fréquence pourrait être imposée, et un bon candidat semblait être l'équation de crête des valeurs extrêmes "a * exp (-1.0 * exp (-1.0 * ((xb) / c )) - ((xb) / c) + 1.0) + offset "comme indiqué ci-dessous.

Voici un ajusteur de courbe Python graphique pour cette équation, en haut du fichier, je charge les données que j'ai extraites afin que vous besoin de le remplacer par les données réelles. Cet ajusteur utilise le module d'algorithme génétique différentielle_évolution de scipy pour estimer les valeurs initiales des paramètres de l'ajusteur non linéaire, qui utilise l'algorithme Latin Hypercube pour assurer une recherche approfondie de l'espace des paramètres et nécessite des limites dans lesquelles la recherche. Ici, ces limites sont prises à partir des valeurs maximum et minimum des données.

En soustrayant les prédictions du modèle de cette courbe ajustée, vous ne devriez avoir que la composante sinusoïdale à modéliser. J'ai noté qu'il semble y avoir un pic supplémentaire étroit de faible amplitude à environ x = 275.

amplitude = -1.0362957093184177E+00
center = 3.6632754608370377E+01
width = 5.0813421718648293E+00
Offset = 5.1940843481496088E+00

pi = 3.14159265358979323846 # constant not fitted

y = amplitude * sin(pi * (x - center) / width) + Offset

MISE À JOUR ------- Si la composante sinusoïdale haute fréquence est constante (ce que je ne sais pas), la modélisation d'une petite partie des données avec seulement quelques cycles sera suffisante pour déterminer l'équation et les estimations des paramètres initiaux pour ajuster la partie d'onde sinusoïdale du modèle. Ici, j'ai fait cela avec le résultat suivant:

à partir de l'équation suivante:

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)

La combinaison de ces deux modèles en utilisant les données réelles, plutôt que mes données extraites par nuage de points, devrait être proche de ce dont vous avez besoin.

J'ai enfin trouvé un moyen de modéliser les ondulations et je peux répondre à ma propre question. Cet 2006 article fait des ajustements de courbes sur des ondulations qui ressemblent à mon ensemble de données.

Tout d'abord, j'ai fait un ajustement polynomial des moindres carrés, puis j'ai soustrait cette courbe polynomiale des données d'origine. Cela ne m'a laissé que les ondulations. En appliquant la transformée de Fourier, j'ai choisi les fréquences dominantes qui me permettent de reconstruire les ondulations sinusoïdales. Ensuite, j'ai simplement ajouté ces ondulations à la courbe polynomiale que j'avais obtenue au début. Cela l'a fait.