Curva di decadimento esponenziale in numpy e scipy

Question

Sto avendo un po 'di problemi con l'adattamento di una curva ad alcuni dati, ma non riesco a capire dove sto andando male.

In passato ho fatto questo con numpy.linalg.lstsq per le funzioni esponenziali e scipy.optimize.curve_fit per le funzioni sigma. Questa volta ho voluto creare uno script che mi permettesse di specificare varie funzioni, determinare i parametri e testarne la compatibilità con i dati. Mentre facevo questo, ho notato che Scipy e Numpy lstsq sembrano fornire risposte diverse per lo stesso insieme di dati e la stessa funzione. La funzione è semplicemente y = e^(l*x) ed è vincolata in modo tale che y=1 a x=0.

La linea di tendenza di Excel è d'accordo con il risultato di Numpy lstsq, ma siccome Scipy leastsq è in grado di svolgere qualsiasi funzione, sarebbe opportuno capire qual è il problema.

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

## Sampled data 
x = np.array([0, 14, 37, 975, 2013, 2095, 2147]) 
y = np.array([1.0, 0.764317544, 0.647136491, 0.070803763, 0.003630962,  0.001485394,  0.000495131]) 

# function 
fp = lambda p, x: np.exp(p*x) 

# error function 
e = lambda p, x, y: (fp(p, x) - y) 

# using scipy least squares 
l1, s = optimize.leastsq(e, -0.004, args=(x,y)) 
print l1 
# [-0.0132281] 


# using numpy least squares 
l2 = np.linalg.lstsq(np.vstack([x, np.zeros(len(x))]).T,np.log(y))[0][0] 
print l2 
# -0.00313461628963 (same answer as Excel trend line) 

# smooth x for plotting 
x_ = np.arange(0, x[-1], 0.2) 

plt.figure() 
plt.plot(x, y, 'rx', x_, fp(l1, x_), 'b-', x_, fp(l2, x_), 'g-') 
plt.show() 

Modifica - ulteriori informazioni

La MWE sopra comprende un piccolo campione del set di dati. Quando si adattano i dati effettivi la curva scipy.optimize.curve_fit presenta un R^2 di 0,82, mentre la curva numpy.linalg.lstsq, che è la stessa calcolata da Excel, ha un R^2 di 0,41 .

Answer 1

Si stanno riducendo al minimo le diverse funzioni di errore.

Quando si utilizza numpy.linalg.lstsq, la funzione di errore sia minimizzato è

np.sum((np.log(y) - p * x)**2) 

mentre scipy.optimize.leastsq minimizza la funzione

np.sum((y - np.exp(p * x))**2) 

Il primo caso richiede una dipendenza lineare tra le variabili dipendenti ed indipendenti, ma la la soluzione è conosciuta analiticamente, mentre la seconda può gestire qualsiasi dipendenza, ma si basa su un metodo iterativo.

In una nota separata, non posso provarlo in questo momento, ma quando si utilizza numpy.linalg.lstsq, ti non c'è bisogno di vstack una fila di zeri, le seguenti opere così:

l2 = np.linalg.lstsq(x[:, None], np.log(y))[0][0] 

Answer 2

Per esporre un po 'sul punto di Jaime, qualsiasi trasformazione non lineare dei dati porterà a una funzione di errore diversa e quindi a soluzioni diverse. Ciò porterà a intervalli di confidenza diversi per i parametri di adattamento. Quindi hai tre possibili criteri da utilizzare per prendere una decisione: quale errore vuoi minimizzare, quali parametri vuoi più sicuro, e infine, se stai usando il fitting per prevedere un certo valore, quale metodo produce meno errori nell'intervallo valore previsto Giocare un po 'analiticamente e in Excel suggerisce che diversi tipi di rumore nei dati (ad esempio se la funzione di rumore scala l'ampiezza, influenza la costante di tempo o è additivo) porta a diverse scelte di soluzione.

Aggiungo anche che mentre questo trucco "funziona" per il decadimento esponenziale a 0, non può essere utilizzato nel caso più generale (e comune) di esponenziali smorzati (in aumento o in calo) a valori che non possono essere supposto essere 0.