L'ajustement d'une courbe à une distribution en loi de puissance avec curve_fit ne fonctionne pas

Question

J'essaie de trouver une courbe qui convienne à mes données qui, visuellement, semblent avoir une distribution de type loi de puissance.

J'espérais pouvoir utiliser scipy.optimize.curve_fit, mais quelle que soit la fonction ou la normalisation des données que j'essaie, j'obtiens soit une RuntimeError (paramètres non trouvés ou dépassement de capacité), soit une courbe qui ne correspond pas du tout à mes données. Veuillez m'aider à comprendre ce que je fais de mal ici.

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

df = pd.DataFrame({
            'x': [ 1000, 3250, 5500, 10000, 32500, 55000, 77500, 100000, 200000 ],
            'y': [ 1100, 500, 288, 200, 113, 67, 52, 44, 5 ]
        })
df.plot(x='x', y='y', kind='line', style='--ro', figsize=(10, 5))

def func_powerlaw(x, m, c, c0):
    return c0 + x**m * c

target_func = func_powerlaw

X = df['x']
y = df['y']

popt, pcov = curve_fit(target_func, X, y)

plt.figure(figsize=(10, 5))
plt.plot(X, target_func(X, *popt), '--')
plt.plot(X, y, 'ro')
plt.legend()
plt.show()

Sortie

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-243-17421b6b0c14> in <module>()
     18 y = df['y']
     19 
---> 20 popt, pcov = curve_fit(target_func, X, y)
     21 
     22 plt.figure(figsize=(10, 5))

/Users/evgenyp/.virtualenvs/kindle-dev/lib/python2.7/site-packages/scipy/optimize/minpack.pyc in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, **kwargs)
    653         cost = np.sum(infodict['fvec'] ** 2)
    654         if ier not in [1, 2, 3, 4]:
--> 655             raise RuntimeError("Optimal parameters not found: " + errmsg)
    656     else:
    657         res = least_squares(func, p0, args=args, bounds=bounds, method=method,

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 800.

Answer 1

Comme l'indique l'historique, le nombre maximal d'évaluations de la fonction a été atteint sans trouver de point stationnaire (pour mettre fin à l'algorithme). Vous pouvez augmenter le nombre maximum en utilisant l'option maxfev . Pour cet exemple, le réglage maxfev=2000 est suffisamment grande pour que l'algorithme se termine avec succès.

Cependant, la solution n'est pas satisfaisante. Ceci est dû au fait que l'algorithme choisit une estimation initiale (par défaut) des variables, qui, pour cet exemple, n'est pas bonne (le grand nombre d'itérations nécessaires en est un indicateur). En fournissant un autre point d'initialisation (trouvé par simple essai et erreur), on obtient un bon ajustement, sans qu'il soit nécessaire d'augmenter le nombre d'itérations. maxfev .

Les deux ajustements et une comparaison visuelle avec les données sont présentés ci-dessous.

x = np.asarray([ 1000, 3250, 5500, 10000, 32500, 55000, 77500, 100000, 200000 ])
y = np.asarray([ 1100, 500, 288, 200, 113, 67, 52, 44, 5 ])

sol1 = curve_fit(func_powerlaw, x, y, maxfev=2000 )
sol2 = curve_fit(func_powerlaw, x, y, p0 = np.asarray([-1,10**5,0]))