Come forzare l'intercettazione zero nella regressione lineare?

Question

Sono un po 'un novello, quindi mi scuso se questa domanda ha già avuto risposta, ho dato un'occhiata e non ho potuto trovare esattamente quello che stavo cercando.

Ho alcuni dati più o meno lineari della forma

x = [0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 20.0, 40.0, 60.0, 80.0] 
y = [0.50505332505407008, 1.1207373784533172, 2.1981844719020001, 3.1746209003398689, 4.2905482471260044, 6.2816226678076958, 11.073788414382639, 23.248479770546009, 32.120462301367183, 44.036117671229206, 54.009003143831116, 102.7077685684846, 185.72880217806673, 256.12183145545811, 301.97120103079675] 

Sto usando scipy.optimize.leastsq per adattare una regressione lineare a questo:

def lin_fit(x, y): 
    '''Fits a linear fit of the form mx+b to the data''' 
    fitfunc = lambda params, x: params[0] * x + params[1] #create fitting function of form mx+b 
    errfunc = lambda p, x, y: fitfunc(p, x) - y    #create error function for least squares fit 

    init_a = 0.5       #find initial value for a (gradient) 
    init_b = min(y)       #find initial value for b (y axis intersection) 
    init_p = numpy.array((init_a, init_b)) #bundle initial values in initial parameters 

    #calculate best fitting parameters (i.e. m and b) using the error function 
    p1, success = scipy.optimize.leastsq(errfunc, init_p.copy(), args = (x, y)) 
    f = fitfunc(p1, x)   #create a fit with those parameters 
    return p1, f  

e funziona in modo bello (anche se non sono certo se scipy.optimize è la cosa giusta da usare qui, potrebbe essere un po 'esagerato?).

Tuttavia, a causa del modo in cui i punti di dati si trovano, non mi dà un'intercettazione sull'asse y a 0. So però che deve essere zero in questo caso, if x = 0 than y = 0.

Esiste un modo per forzarlo?

Answer 1

Non sono esperto in questi moduli, ma ho una certa esperienza in statistica, quindi ecco quello che vedo. È necessario modificare la funzione in forma da

fitfunc = lambda params, x: params[0] * x + params[1] 

a:

fitfunc = lambda params, x: params[0] * x 

rimuovere Anche la linea:

init_b = min(y) 

e cambiare la riga successiva a:

init_p = numpy.array((init_a)) 

Questo dovrebbe sbarazzarsi del secondo parametro questo sta producendo l'intercetta y e passa la linea adattata attraverso l'origine. Potrebbero esserci altre modifiche minori che potresti dover fare per questo nel resto del tuo codice.

Ma sì, non sono sicuro che questo modulo funzionerà se si rimuove il secondo parametro in questo modo. Dipende dal funzionamento interno del modulo se può accettare questa modifica. Ad esempio, non so dove sia stato inizializzato lo params, l'elenco di parametri, quindi non so se fare proprio questo cambierà la sua lunghezza.

E per inciso, dato che hai menzionato, questo in realtà penso sia un modo un po 'esagerato per ottimizzare solo una pendenza. È possibile leggere la regressione lineare un po 'e scrivere un piccolo codice per farlo da soli dopo un back-of-the-envelope. È abbastanza semplice e diretto, davvero. In effetti, ho appena fatto alcuni calcoli, e immagino che la pendenza ottimizzata sarà solo <xy>/<x^2>, cioè la media dei prodotti x * y divisa per la media di x^2.

Answer 2

Come @AbhranilDas menzionato, basta usare un metodo lineare. Non è necessario un risolutore non lineare come scipy.optimize.lstsq.

In genere, utilizza numpy.polyfit per adattare una linea per i dati, ma in questo caso è necessario fare uso numpy.linalg.lstsq direttamente, come si desidera impostare l'intercetta a zero.

Come un esempio veloce:

import numpy as np 
import matplotlib.pyplot as plt 

x = np.array([0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 
       20.0, 40.0, 60.0, 80.0]) 

y = np.array([0.50505332505407008, 1.1207373784533172, 2.1981844719020001, 
       3.1746209003398689, 4.2905482471260044, 6.2816226678076958, 
       11.073788414382639, 23.248479770546009, 32.120462301367183, 
       44.036117671229206, 54.009003143831116, 102.7077685684846, 
       185.72880217806673, 256.12183145545811, 301.97120103079675]) 

# Our model is y = a * x, so things are quite simple, in this case... 
# x needs to be a column vector instead of a 1D vector for this, however. 
x = x[:,np.newaxis] 
a, _, _, _ = np.linalg.lstsq(x, y) 

plt.plot(x, y, 'bo') 
plt.plot(x, a*x, 'r-') 
plt.show()