Scipy.optimize.curve_fit ne correspondra pas à la loi de puissance du cosinus

Question

Depuis plusieurs heures maintenant, j'essaie d'adapter un modèle à un ensemble de données (généré) en tant que casus pour un problème avec lequel j'ai lutté. J'ai généré des points de données pour la fonction f (x) = A * cos^n (x) + b, et j'ai ajouté du bruit. Lorsque je tente d'adapter l'ensemble de données avec cette fonction et curve_fit, je reçois l'erreur

./tester.py:10: RuntimeWarning: invalid value encountered in power 
return Amp*(np.cos(x))**n + b 
/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py:690: OptimizeWarning: Covariance of the parameters could not be estimated category=OptimizeWarning) 

Le code que je utilise pour générer les points de données et adapter le modèle est le suivant:

#!/usr/bin/env python 

from __future__ import print_function 
import numpy as np 
from scipy.optimize import curve_fit 
from matplotlib.pyplot import figure, show, rc, plot 

def f(x, Amp, n, b): 

    return np.real(Amp*(np.cos(x))**n + b) 

x = np.arange(0, 6.28, 0.01) 
randomPart = np.random.rand(len(x))-0.5 
fig = figure() 
sample = f(x, 5, 2, 5)+randomPart 
frame = fig.add_subplot(1,1,1) 

frame.plot(x, sample, label="Sample measurements") 

popt, pcov = curve_fit(f, x, sample, p0=(1,1,1)) 

modeldata = f(x, popt[0], popt[1], popt[2]) 
print(modeldata) 
frame.plot(x, modeldata, label="Best fit") 

frame.legend() 
frame.set_xlabel("x") 
frame.set_ylabel("y") 

show() 

Le Les données bruyantes sont affichées - voir l'image ci-dessous.

Est-ce que quelqu'un d'entre vous la moindre idée de ce qui se passe? Je soupçonne qu'il a quelque chose à voir avec la loi de puissance entrant dans le domaine complexe, comme la partie réelle de la fonction est nowhere divergent. J'ai essayé de retourner seulement la partie réelle de la fonction, en définissant des limites réalistes dans curve_fit et en utilisant déjà un tableau numpy au lieu d'une liste python pour p0. Je cours la dernière version de scipy disponible, scipy 0.17.0-1.

Answer 1

Le problème est le suivant:

>>> (-2)**1.1 
(-2.0386342710747223-0.6623924280875919j) 
>>> np.array(-2)**1.1 
__main__:1: RuntimeWarning: invalid value encountered in power 
nan 

Contrairement à flotteurs python indigènes, numpy double refusent généralement de prendre part à des opérations conduisant à des résultats complexes:

>>> np.sqrt(-1) 
__main__:1: RuntimeWarning: invalid value encountered in sqrt 
nan 

En tant que solution rapide, je suggère d'ajouter un appel np.abs à votre fonction, et en utilisant les limites appropriées pour l'ajustement pour s'assurer que cela ne donne pas un faux ajustement. Si votre modèle est proche de la vérité et que votre échantillon (je veux dire le cosinus de votre échantillon) est positif, alors ajouter une valeur absolue autour de ce modèle devrait être non-op (mise à jour: je réalise que ce n'est jamais le cas au dessous de).

def f(x, Amp, n, b): 

    return Amp*(np.abs(np.cos(x)))**n + b # only change here 

Avec ce petit changement que j'obtiens ceci:

Pour référence, les paramètres de l'ajustement sont (4.96482314, 2.03690954, 5.03709923]) comparant à la génération avec (5,2,5).

---

Après avoir un peu plus pensé que je compris que le cosinus sera toujours être négatif pour la moitié de votre domaine (duh). Donc la solution de contournement que j'ai suggérée pourrait être un peu problématique, ou au moins sa correction n'est pas triviale. D'autre part, en pensant à votre formule originale contenant cos(x)^n, avec des valeurs négatives pour cos(x) cela n'a de sens en tant que modèle si n est un nombre entier, sinon vous obtenez un résultat complexe. Puisque nous ne pouvons pas résoudre Diophantine problèmes d'ajustement, nous devons gérer cela correctement.La manière la plus appropriée (par laquelle je veux dire la façon la moins susceptible de biaiser vos données) est la suivante: d'abord faire l'ajustement avec un modèle qui convertit vos données en nombres complexes puis prend la magnitude complexe sur la sortie:

def f(x, Amp, n, b): 

    return Amp*np.abs(np.cos(x.astype(np.complex128))**n) + b 

Ceci est évidemment beaucoup moins efficace que ma solution de contournement, puisque dans chaque étape de montage, nous créons un nouveau maillage, et de faire un travail supplémentaire à la fois sous la forme d'arithmétique complexe et un calcul de magnitude supplémentaire. Cela me donne l'ajustement suivant, même sans limites définies:

Les paramètres sont (5.02849409, 1.97655728, 4.96529108). Ceux-ci sont proches aussi. Cependant, si nous remettons ces valeurs dans le modèle réel (sans np.abs), nous obtenons des parties imaginaires aussi grandes que -0.37, ce qui n'est pas écrasant mais significatif. Donc, la deuxième étape devrait être de refaire l'ajustement avec un bon modèle --- celui qui a un exposant entier. Prenez l'exposant 2 qui est évident de votre ajustement, et faire un nouvel ajustement avec ce modèle. Je ne crois pas que toute autre approche vous donne un résultat mathématiquement solide. Vous pouvez également commencer à partir de l'original popt, en espérant que c'est en effet proche de la vérité. Bien sûr, nous pourrions utiliser la fonction d'origine avec certains curry, mais il est beaucoup plus rapide d'utiliser une version spécifique double spécifique de votre modèle. Notez que j'ai modifié quelques éléments dans votre code, ce qui n'affecte pas vraiment mon point de vue. La figure de ce qui précède, de sorte que celui qui montre à la fois l'ajustement auxiliaire et le bon un:

La sortie:

Auxiliary fit parameters: [ 5.02628994 2.00886409 5.00652371] 
Final fit parameters: [5.0288141074549699, 2.0, 5.0009730316739462] 

simplement répéter: alors qu'il pourrait y avoir aucune différence visuelle entre l'ajustement auxiliaire et le bon, seul ce dernier donne une réponse significative à votre problème.