J'ai une routine de référence Matlab que j'essaie de convertir en numpy/scipy. J'ai rencontré un problème d'ajustement de courbe que je ne peux pas résoudre en Python. Voici donc un exemple simple qui illustre le problème. Les données sont complètement synthétiques et ne font pas partie du problème.

Disons que j'essaie d'adapter un modèle linéaire de données bruitées -

x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
y = [0.1075, 1.3668, 1.5482, 3.1724, 4.0638, 4.7385, 5.9133, 7.0685, 8.7157, 9.5539]

Pour la solution non pondérée dans Matlab, je coderais

g = @(m, b, x)(m*x + b)
f = fittype(g)
bestfit = fit(x, y, g)

Ce qui produit une solution de bestfit.m = 1.048, bestfit.b = -0.09219

L'exécution de ces données via scipy.optimize.curve_fit() produit des résultats identiques.

Si à la place, l'ajustement utilise une fonction de décroissance pour réduire l'impact des points de données

dw = [0.7290, 0.5120, 0.3430, 0.2160, 0.1250, 0.0640, 0.0270, 0.0080,  0.0010, 0]
weightedfit = fit(x, y, g, 'Weights', dw)

Cela produit une pente si 0,944 et un décalage de 0,1484.

Je n'ai pas compris comment évoquer ce résultat à partir de scipy.optimize.curve_fit en utilisant le paramètre sigma. Si je transmets les poids fournis à Matlab, le '0' provoque une exception de division par zéro. Il est clair que Matlab et scipy pensent très différemment à la signification des poids dans la routine d'optimisation sous-jacente. Existe-t-il un moyen simple de convertir entre les deux qui me permette de fournir une fonction de pondération qui produit des résultats identiques ?

0
bbeauchaine87408 21 nov. 2019 à 23:20

1 réponse

Meilleure réponse

Ok, donc après une enquête plus approfondie, je peux offrir la réponse, au moins pour cet exemple simple.

import numpy as np
import scipy as sp
import scipy.optimize

def modelFun(x, m, b):
    return m * x + b

def testFit():
    w = np.diag([1.0, 1/0.7290, 1/0.5120, 1/0.3430, 1/0.2160, 1/0.1250, 1/0.0640, 1/0.0270, 1/0.0080, 1/0.0010])
    x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    y = np.array([0.1075, 1.3668, 1.5482, 3.1724, 4.0638, 4.7385, 5.9133, 7.0685, 8.7157, 9.5539])

    popt = sp.optimize.curve_fit(modelFun, x, y, sigma=w)

    print(popt[0])
    print(popt[1])

Ce qui produit le résultat souhaité.

Afin de forcer sp.optimize.curve_fit à minimiser la même métrique chisq que Matlab à l'aide de la boîte à outils d'ajustement de courbe, vous devez faire deux choses :

  1. Utiliser l'inverse des facteurs de poids
  2. Créez une matrice diagonale à partir des nouveaux facteurs de pondération. Selon la référence scipy :

sigma Aucun ou séquence de longueur M ou tableau MxM, facultatif Détermine l'incertitude dans les données y. Si nous définissons les résidus comme r = ydata - f(xdata, *popt), alors l'interprétation de sigma dépend de son nombre de dimensions :

Un sigma 1-d doit contenir des valeurs d'écarts types d'erreurs dans les données y. Dans ce cas, la fonction optimisée est chisq = sum((r / sigma) ** 2).

Un sigma 2D doit contenir la matrice de covariance des erreurs dans les données y. Dans ce cas, la fonction optimisée est chisq = r.T @ inv(sigma) @ r.

Nouveau dans la version 0.19.

Aucun (par défaut) équivaut à un sigma 1-D rempli de uns.

1
bbeauchaine87408 3 déc. 2019 à 15:41