Perché scipy.optimize.curve_fit non si adatta ai dati?

Question

Sto provando ad adattare un esponenziale ad alcuni dati per un po 'usando scipy.optimize.curve_fit ma sto avendo delle vere difficoltà. Non riesco davvero a vedere alcun motivo per cui questo non funzionerebbe, ma produce solo una linea di confine, nessuna idea del perché!

Tutto l'aiuto sarebbe apprezzato

from __future__ import division 
import numpy 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as pyplot 

def func(x,a,b,c): 
    return a*numpy.exp(-b*x)-c 


yData = numpy.load('yData.npy') 
xData = numpy.load('xData.npy') 

trialX = numpy.linspace(xData[0],xData[-1],1000) 

# Fit a polynomial 
fitted = numpy.polyfit(xData, yData, 10)[::-1] 
y = numpy.zeros(len(trailX)) 
for i in range(len(fitted)): 
    y += fitted[i]*trialX**i 

# Fit an exponential 
popt, pcov = curve_fit(func, xData, yData) 
yEXP = func(trialX, *popt) 

pyplot.figure() 
pyplot.plot(xData, yData, label='Data', marker='o') 
pyplot.plot(trialX, yEXP, 'r-',ls='--', label="Exp Fit") 
pyplot.plot(trialX, y, label = '10 Deg Poly') 
pyplot.legend() 
pyplot.show() 

xData = [1e-06, 2e-06, 3e-06, 4e-06, 
5e-06, 6e-06, 7e-06, 8e-06, 
9e-06, 1e-05, 2e-05, 3e-05, 
4e-05, 5e-05, 6e-05, 7e-05, 
8e-05, 9e-05, 0.0001, 0.0002, 
0.0003, 0.0004, 0.0005, 0.0006, 
0.0007, 0.0008, 0.0009, 0.001, 
0.002, 0.003, 0.004, 0.005, 
0.006, 0.007, 0.008, 0.009, 0.01] 

yData = 
[6.37420666067e-09, 1.13082012115e-08, 
1.52835756975e-08, 2.19214493931e-08, 2.71258852882e-08, 3.38556130078e-08, 3.55765277358e-08, 
4.13818145846e-08, 4.72543475372e-08, 4.85834751151e-08, 9.53876562077e-08, 1.45110636413e-07, 
1.83066627931e-07, 2.10138415308e-07, 2.43503982686e-07, 2.72107045549e-07, 3.02911771395e-07, 
3.26499455951e-07, 3.48319349445e-07, 5.13187669283e-07, 5.98480176303e-07, 6.57028222701e-07, 
6.98347073045e-07, 7.28699930335e-07, 7.50686502279e-07, 7.7015576866e-07, 7.87147246927e-07, 
7.99607141001e-07, 8.61398763228e-07, 8.84272900407e-07, 8.96463883243e-07, 9.04105135329e-07, 
9.08443443149e-07, 9.12391264185e-07, 9.150842683e-07, 9.16878548643e-07, 9.18389990067e-07] 

Answer 1

algoritmi numerici tendono a lavorare meglio quando non sono alimentati estremamente piccole (o grandi) numeri.

In questo caso, il grafico mostra i valori x e y estremamente piccoli. Se li si scala, l'adattamento è notevole meglio:

xData = np.load('xData.npy')*10**5 
yData = np.load('yData.npy')*10**5 

---

from __future__ import division 

import os 
os.chdir(os.path.expanduser('~/tmp')) 

import numpy as np 
import scipy.optimize as optimize 
import matplotlib.pyplot as plt 

def func(x,a,b,c): 
    return a*np.exp(-b*x)-c 


xData = np.load('xData.npy')*10**5 
yData = np.load('yData.npy')*10**5 

print(xData.min(), xData.max()) 
print(yData.min(), yData.max()) 

trialX = np.linspace(xData[0], xData[-1], 1000) 

# Fit a polynomial 
fitted = np.polyfit(xData, yData, 10)[::-1] 
y = np.zeros(len(trialX)) 
for i in range(len(fitted)): 
    y += fitted[i]*trialX**i 

# Fit an exponential 
popt, pcov = optimize.curve_fit(func, xData, yData) 
print(popt) 
yEXP = func(trialX, *popt) 

plt.figure() 
plt.plot(xData, yData, label='Data', marker='o') 
plt.plot(trialX, yEXP, 'r-',ls='--', label="Exp Fit") 
plt.plot(trialX, y, label = '10 Deg Poly') 
plt.legend() 
plt.show() 

Nota che dopo ridimensionamento xData e yData, i parametri restituiti da curve_fit devono essere riscalati. In questo caso, a, b e c ciascuno devono essere divisi per 10 ** 5 per ottenere i parametri corretti per i dati originali.

---

Un'obiezione che potresti avere sopra è che il ridimensionamento deve essere scelto piuttosto "attentamente". (Leggi: non tutte le scelte ragionevoli di scala funziona!)

È possibile migliorare la robustezza di curve_fit fornendo un'ipotesi iniziale ragionevole per i parametri. Di solito hai qualche a priori conoscenza dei dati che possono motivare ipotesi di tipo "back-of-the-envelope" per valori di parametro ragionevoli.

Ad esempio, la chiamata curve_fit con

guess = (-1, 0.1, 0) 
popt, pcov = optimize.curve_fit(func, xData, yData, guess) 

aiuta a migliorare la gamma di scale su cui curve_fit riesce in questo caso.

Answer 2

Un (lieve) miglioramento di questa soluzione, non tenendo conto della conoscenza a priori dei dati potrebbe essere il seguente: Prendere la media inversa del set di dati e utilizzarla come "fattore di scala" da passare al leastsq() sottostante chiamato da curve_fit(). Ciò consente all'installatore di funzionare e restituisce i parametri sulla scala originale dei dati.

La linea di riferimento è:

popt, pcov = curve_fit(func, xData, yData) 

che diventa:

popt, pcov = curve_fit(func, xData, yData, 
    diag=(1./xData.mean(),1./yData.mean())) 

Ecco l'esempio completo che produce questa immagine:

from __future__ import division 
import numpy 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as pyplot 

def func(x,a,b,c): 
    return a*numpy.exp(-b*x)-c 


xData = numpy.array([1e-06, 2e-06, 3e-06, 4e-06, 5e-06, 6e-06, 
7e-06, 8e-06, 9e-06, 1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 
7e-05, 8e-05, 9e-05, 0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 
0.0006, 0.0007, 0.0008, 0.0009, 0.001, 0.002, 0.003, 0.004, 0.005 
, 0.006, 0.007, 0.008, 0.009, 0.01]) 

yData = numpy.array([6.37420666067e-09, 1.13082012115e-08, 
1.52835756975e-08, 2.19214493931e-08, 2.71258852882e-08, 
3.38556130078e-08, 3.55765277358e-08, 4.13818145846e-08, 
4.72543475372e-08, 4.85834751151e-08, 9.53876562077e-08, 
1.45110636413e-07, 1.83066627931e-07, 2.10138415308e-07, 
2.43503982686e-07, 2.72107045549e-07, 3.02911771395e-07, 
3.26499455951e-07, 3.48319349445e-07, 5.13187669283e-07, 
5.98480176303e-07, 6.57028222701e-07, 6.98347073045e-07, 
7.28699930335e-07, 7.50686502279e-07, 7.7015576866e-07, 
7.87147246927e-07, 7.99607141001e-07, 8.61398763228e-07, 
8.84272900407e-07, 8.96463883243e-07, 9.04105135329e-07, 
9.08443443149e-07, 9.12391264185e-07, 9.150842683e-07, 
9.16878548643e-07, 9.18389990067e-07]) 

trialX = numpy.linspace(xData[0],xData[-1],1000) 

# Fit a polynomial 
fitted = numpy.polyfit(xData, yData, 10)[::-1] 
y = numpy.zeros(len(trialX)) 
for i in range(len(fitted)): 
    y += fitted[i]*trialX**i 

# Fit an exponential 
popt, pcov = curve_fit(func, xData, yData, 
    diag=(1./xData.mean(),1./yData.mean())) 
yEXP = func(trialX, *popt) 

pyplot.figure() 
pyplot.plot(xData, yData, label='Data', marker='o') 
pyplot.plot(trialX, yEXP, 'r-',ls='--', label="Exp Fit") 
pyplot.plot(trialX, y, label = '10 Deg Poly') 
pyplot.legend() 
pyplot.show() 

Answer 3

modello a*exp(-b*x)+c adatta bene i dati, ma vi suggerisco una piccola modifica:
uso questo, invece

a*x*exp(-b*x)+c

buona fortuna