Método de mínimos cuadrados
Veremos a continuación como emplear el método de mínimos cuadrados para realizar el ajuste de un modelo logístico a valores de concentración de bacterias medidos experimentalmente.
Datos disponibles
Se ha realizado una serie de experimentos en los que se ha evaluado la concentración de bacterias (\(y\)) en función del tiempo transcurrido (\(x\)) medido en minutos. Los datos disponibles son:
x | y |
0 | 0.068 |
5 | 0.127 |
10 | 0.077 |
15 | 0.359 |
20 | 0.738 |
25 | 0.885 |
30 | 0.828 |
35 | 1.058 |
40 | 0.989 |
45 | 1.122 |
50 | 1.090 |
Representamos estos datos gráficamente:
Modelo logístico para la concentración bacteriana.
La función logística a la que se pretende ajustar estos datos es de la forma:
\[f(x)=\frac{e^{b_1+b_2x}}{1+e^{b_1+b_2x}}\]
donde los valores de \(b_1\) y \(b_2\) no se conocen.
Método de mínimos cuadrados
Para obtener los valores de los parámetros \(b_1\) y \(b_2\) aplicaremos el método de mínimos cuadrados consistente en calcular \(b_1\) y \(b_2\) como los valores que minimizan la suma de los cuadrados de las distancias entre los valores observados \(y_i\) y los valores predichos \(\hat{y}_i\) por el modelo:
\[\min \sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2}\]
Para cada \(y_i\) el valor predicho se obtiene como: \[\hat{y}_{i}=f\left(x_i\right)=\frac{e^{b_1+b_2x_{i}}}{1+e^{b_1+b_2x_{i}}}\] Por tanto, los valores de \(b_1\) y \(b_2\) se obtienen minimizando la función: \[L\left(b_1,b_2\right)=\sum_{i=1}^{n}\left(y_{i}-\frac{e^{b_1+b_2x_{i}}}{1+e^{b_1+b_2x_{i}}}\right)^{2}\] Nótese que en esta función los valores de las \(x_i\) y las \(y_i\) son conocidos y fijos (son los valores obtenidos experimentalmente que figuran en la tabla inicial).
Obtención de los valores de los parámetros mediante el método de mínimos cuadrados.
Para obtener los valores de \(b_1\) y \(b_2\) que minimizan esta función debemos derivar respecto a cada uno de los parámetros e igualar a 0:
\[\begin{cases} \frac{\partial L}{\partial b_1} & =0\\ \frac{\partial L}{\partial b_2} & =0 \end{cases}\]
Tenemos:
\[\frac{\partial L}{\partial b_1} =\sum_{i=1}^{n}2\left(y_{i}-\frac{e^{b_1+b_2x_{i}}}{1+e^{b_1+b_2x_{i}}}\right)\left(-\frac{\left(e^{b_1+b_2x_{i}}\right)\left(1+e^{b_1+b_2x_{i}}\right)-e^{b_1+b_2x_{i}}e^{b_1+b_2x_{i}}}{\left(1+e^{b_1+b_2x_{i}}\right)^{2}}\right)=\] \[=-2\sum_{i=1}^{n}\left(y_{i}-\frac{e^{b_1+b_2x_{i}}}{1+e^{b_1+b_2x_{i}}}\right)\left(\frac{\left(e^{b_1+b_2x_{i}}\right)}{\left(1+e^{b_1+b_2x_{i}}\right)^{2}}\right)=-2\sum_{i=1}^{n}\left(y_{i}-f\left(x_{i}\right)\right)\frac{f\left(x_{i}\right)}{1+e^{b_1+b_2x_{i}}}\]
Asimismo:
\[\frac{\partial L}{\partial b_2} =\sum_{i=1}^{n}2\left(y_{i}-\frac{e^{b_1+b_2x_{i}}}{1+e^{b_1+b_2x_{i}}}\right)\left(-\frac{x_{i}\left(e^{b_1+b_2x_{i}}\right)\left(1+e^{b_1+b_2x_{i}}\right)-e^{b_1+b_2x_{i}}x_{i}e^{b_1+b_2x_{i}}}{\left(1+e^{b_1+b_2x_{i}}\right)^{2}}\right)=\] \[=-2\sum_{i=1}^{n}\left(y_{i}-\frac{e^{b_1+b_2x_{i}}}{1+e^{b_1+b_2x_{i}}}\right)\left(\frac{x_{i}e^{b_1+b_2x_{i}}}{\left(1+e^{b_1+b_2x_{i}}\right)^{2}}\right)=-2\sum_{i=1}^{n}\left(y_{i}-f\left(x_{i}\right)\right)\frac{x_{i}f\left(x_{i}\right)}{1+e^{b_1+b_2x_{i}}}\]
Por tanto, el sistema de ecuaciones a resolver es:
\[\begin{cases} \sum_{i=1}^{n}\left(y_{i}-f\left(x_{i}\right)\right)\frac{f\left(x_{i}\right)}{1+e^{b_1+b_2x_{i}}} & =0\\ \sum_{i=1}^{n}\left(y_{i}-f\left(x_{i}\right)\right)\frac{x_{i}f\left(x_{i}\right)}{1+e^{b_1+b_2x_{i}}} & =0 \end{cases}\]
Llamando:
\[M_1\left(b_1,b_2\right)=\sum_{i=1}^{n}\left(y_{i}-f\left(x_{i}\right)\right)\frac{f\left(x_{i}\right)}{1+e^{b_1+b_2x_{i}}}\]
\[M_2\left(b_1,b_2\right)=\sum_{i=1}^{n}\left(y_{i}-f\left(x_{i}\right)\right)\frac{x_i f\left(x_{i}\right)}{1+e^{b_1+b_2x_{i}}}\]
y
\[M\left(b_{1},b_{2}\right)=\left(M_{1}\left(b_{1},b_{2}\right),M_{2}\left(b_{1},b_{2}\right)\right)\]
nuestro problema se reduce a encontrar el vector \(\left(b_{1},b_{2}\right)\) para el cual:
\[M\left(b_{1},b_{2}\right)=\left(0,0\right)\]
Implementación en matlab
Para obtener los valores de \(b_1\) y \(b_2\) que resuelven el sistema de ecuaciones anterior, utilizaremos la función fsolve de matlab. Para ello en primer lugar debemos implementar la función vectorial \(M\) anterior, que puede hacerse en matlab del siguiente modo:
datos=readtable("xy1.csv")
x=datos.x
y=datos.y
d=@(b) 1+exp(b(1)+b(2)*x)
f=@(b) exp(b(1)+b(2)*x)./d(b)
M1=@(b) sum(((y-f(b)).*f(b))./d(b))
M2=@(b) sum(((y-f(b)).*x.*f(b))./d(b))
M=@(b) [M1(b), M2(b)]
M([1,1])
Obtención de una aproximación inicial a la solución
Para resolver el sistema de ecuaciones, es necesario proporcionar a Matlab una solución inicial. Para obtener valores iniciales de \(b_1\) y \(b_2\) observamos que:
- Para \(x=0\) la \(y\) vale 0.06. Por tanto \(f(0)\simeq 0.07\). Como: \[f\left(0\right)=\frac{e^{b_{0}}}{1+e^{b_{0}}}\]
hacemos:
\[\frac{e^{b_{0}}}{1+e^{b_{0}}}=0.07\]
de donde:
\[0.93e^{b_{0}}\simeq0.07\Rightarrow b_{0}\simeq\log\left(\frac{0.07}{0.93}\right)\simeq-2.6\]
Asimismo, vemos que para valores de \(x\) por encima de 40, los valores de \(y\) observados están en 1 o ligeramente por encima de 1. Como la función \(f(x)\) tal como está definida no puede tomar valores mayores que 1, aproximamos:
\[f\left(50\right)\simeq0.999\]
de donde:
\[\frac{e^{-2.6+b_{1}\cdot50}}{1+e^{-2.6+b_{1}\cdot50}}\simeq0.999\Longrightarrow0.001e^{-2.6+b_{1}\cdot50}\simeq0.999\Longrightarrow\]
\[-2.6+b_{1}\cdot50\simeq\log\left(\frac{0.999}{0.001}\right)\Longrightarrow b_{1}\simeq\frac{\log\left(999\right)+2.6}{50}=0.19\]
Por tanto, los datos sugieren que una buena aproximación inicial para los parámetros del modelo es \(\left(b_0,b_1\right)=\left(-2.6, 0.19\right)\).
Resolución utilizando fsolve()
Podemos introducir esta aproximación inicial en fsolve mediante:
que produce como resultado:
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
sol =
-4.3443 0.25574Por tanto la solución es \(b_0=-4.3443\), \(b_1=0.25574\). Podemos confirmar que efectivamente es solución del sistema de ecuaciones calculando:
Podemos también dibujar la curva obtenida con estos parámetros:
f=@(b,x) exp(b(1)+b(2)*x)./(1+exp(b(1)+b(2)*x))
plot(x,y,"or")
hold on
xx=0:50
yy=f(sol,xx)
plot(xx,yy,"blue")
hold off
Como vemos, con estos parámetros la función se ajusta razonablemente bien a la nube de puntos observada.
Falsas soluciones
Es preciso tener en cuenta que si partimos de otra aproximación inicial, es posible que matlab encuentre otra “solución” distinta que no resuelva nuestro problema. Por ejemplo, eligiendo como solución inicial \(b=[0,1]\):
devuelve como resultado:
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
sol1 =
-2.5682 2.714Evaluamos el valor de la función \(M\) en este punto:
Como vemos, aunque \(M(sol1) \simeq (0,0)\) la solución no es tan buena (tan próxima a 0) como en el caso anterior. Además, si representamos la gráfica de \(f(x)\) usando los valores \(sol1 = \left(b_1,b_2\right)=\) \(=(-2.5682, 2.714)\), obtenemos:
y como puede apreciarse, la función así obtenida no se ajusta bien a la nube de puntos.
Como segundo ejemplo, si elegimos como valor inicial \(b=[-1,0]\) obtenemos
sol2=fsolve(M,[-1,0])
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
sol2 =
-2.6097 -0.9915Nuevamente evaluamos \(M\) en este punto y comprobamos que aunque el valor de \(M\) es próximo a (0,0), no es tan próximo como con la primera solución:
La representación gráfica con estos valores sería:
y como podemos comprobar, la función resultante tampoco se ajusta para nada a la nube de puntos observada.
Por tanto es crucial partir de una buena aproximación inicial a la solución del problema, y comprobar si el resultado satisface nuestras expectativas (produce una función que ses ajuste bien a la nube de puntos)
Mínimos cuadrados utilizando la función lsqcurvefit
Como hemos visto, el método de mínimos cuadrados nos permite encontrar los valores de los parámetros para los que se consigue un buen ajuste entre una nube de puntos experimental y una función teórica que a priori podemos esperar que represente razonablemente bien la trayectoria de dicha nube de puntos (la relación entre la \(x\) y la \(y\)).
No obstante la aplicación que hemos hecho este método nos ha exigido calcular las derivadas de la función de mínimos cuadrados con respecto a cada uno de los parámetros, igualarlas a cero y resolver el sistema de ecuaciones resultante. En muchos casos ésta puede ser una tarea tediosa, y además susceptible de errores difíciles de detectar (al calcular las derivadas es fácil equivocarse en un signo, en un cociente, en un exponente, …). Matlab dispone de la función lsqcurvefit que permite realizar el ajuste por mínimos cuadrados a una función preespecificada, sin necesidad de calcular derivadas, lo que simplifica notablemente el procedimiento.
En realidad el algoritmo detrás de lsqcurvefit emplea técnicas numéricas que aproximan dichas derivadas (si bien es cierto que en algunos casos dicha aproximación no funciona y el algoritmo lsqcurvefit no funciona); no obstante, al igual que en el ejemplo que hemos visto más arriba, la solución depende fuertemente de los valores iniciales que proporcionemos. Si los valores iniciales son inadecuados, el algoritmo nos puede devolver una solución muy alejada de la solución correcta.
Para aplicar lsqcurvefit hay que proporcionar a matlab:
Una implementación de la función que se quiere ajustar; esta función debe depender de dos vectores: el primero con los parámetros de la función (en el ejemplo anterior \(b=\left(b_1,b_2\right)\)) y el segundo con los valores de la variable independiente \(x\)
Una aproximación inicial de los valores de los parámetros.
Los valores observados de la variable independiente \(x\)
Los valores observados de la variable respuesta \(y\)
En nuestro ejemplo, la función que queremos ajustar es:
\[f(x)=\frac{e^{b_1+b_2x}}{1+e^{b_1+b_2x}}\] que implementamos mediante:
Los valores iniciales de los parámetros \((b_1,b_2)\) son \((-2.6,\; 0.19)\).
Los valores de las variables \(x\) e \(y\) los hemos guardado en datos.x y datos.y.
La llamada a la función lsqcurvefit se realizaría entonces del modo siguiente:
que nos devuelve como resultado:
Local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to
its initial value is less than the value of the function tolerance.
ans =
-4.3417 0.25559Como vemos, es una solución casi idéntica a la que ya habíamos obtenido.
Debemos señalar además que el algoritmo lsqcurvefit es muy eficiente, y aunque los valores iniciales no sean demasiado buenos, en muchos casos (no siempre) es capaz de encontrar la solución correcta. Por ejemplo, con fsolve utilizando [-1,1] como valor inicial de \((b_1,b_2)\), ya vimos que se producía una solución incorrecta. Sin embargo, con lsqcurvefit, la solución que se obtiene partiendo de dicha aproximación inicial es la correcta:
>> lsqcurvefit(f , [-1,1], datos.x, datos.y)
Local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to
its initial value is less than the value of the function tolerance.
ans =
-4.3429 0.25566En cualquier caso, debemos insistir en que, a pesar de que lsqcurvefit es más eficiente, en muchas ocasiones si introducimos una mala aproximación inicial, el algoritmo puede no converger (dar la solución), o puede proporcionar una solución incorrecta. Por ejemplo, usando como valor inicial (10,-5):
>> lsqcurvefit(f , [10,-5], datos.x, datos.y)
Local minimum found.
Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
<stopping criteria details>
ans =
-2.6101 -4.9539Si representamos la función con estos parámetros, vemos que el ajuste conseguido es malo:
