Cómo forzar la intercepción cero en la regresión lineal?

Question

Soy un poco nuevo así que me disculpo si esta pregunta ya ha sido respondida, he echado un vistazo y no he podido encontrar específicamente lo que estaba buscando.

Tengo algunos datos más o menos lineales de la 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] 

estoy usando scipy.optimize.leastsq para ajustar una regresión lineal para esto:

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  

y funciona muy bien (aunque no estoy seguro si scipy.optimize es lo correcto para usar aquí, ¿podría ser un poco exagerado?).

Sin embargo, debido a la forma en que se encuentran los datos, no me da una interceptación del eje y en 0. Sin embargo, sé que tiene que ser cero en este caso, if x = 0 than y = 0.

¿Hay alguna manera de forzar esto?

Answer 1

No soy experto en estos módulos, pero tengo cierta experiencia en estadísticas, así que esto es lo que veo. Usted necesita cambiar su función de ajuste de

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

a:

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

También quite la línea:

init_b = min(y) 

Y cambiar la siguiente línea a:

init_p = numpy.array((init_a)) 

Esto debería deshacerse del segundo parámetro que está produciendo la intersección en y y pasa la línea ajustada a través del origen. Es posible que haya un par de modificaciones menores que deba hacer en el resto de su código.

Pero sí, no estoy seguro de si este módulo funcionará si solo arranca el segundo parámetro de esta manera. Depende del funcionamiento interno del módulo en cuanto a si puede aceptar esta modificación. Por ejemplo, no sé dónde se está inicializando params, la lista de parámetros, por lo que no sé si esto cambiará su longitud.

Y, aparte, como mencionaste, esto realmente creo que es una forma excesiva de optimizar solo una pendiente. Puede leer un regresión lineal un poco y escribir un código pequeño para hacerlo usted mismo después de un cálculo al dorso del sobre. Es bastante simple y directo, realmente. De hecho, acabo de hacer algunos cálculos, y supongo que la pendiente optimizada será <xy>/<x^2>, es decir, la media de los productos x * y dividida por la media de x^2.

Answer 2

Como se mencionó en @AbhranilDas, solo use un método lineal. No hay necesidad de un solucionador no lineal como scipy.optimize.lstsq.

Normalmente, utilizaría numpy.polyfit para ajustar una línea a sus datos, pero en este caso necesitará usar numpy.linalg.lstsq directamente, ya que desea configurar el intercepto a cero.

Como un ejemplo rápido:

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()