regresion - ¿Cómo determino los coeficientes para una línea de regresión lineal en MATLAB?

linea de tendencia matlab

Question

Voy a escribir un programa donde la entrada es un conjunto de datos de puntos 2D y la salida es los coeficientes de regresión de la línea de mejor ajuste al minimizar el error MSE mínimo.

Tengo algunos ejemplos de puntos que me gustaría procesar:

  X      Y
1.00    1.00
2.00    2.00
3.00    1.30
4.00    3.75
5.00    2.25

¿Cómo haría esto en MATLAB?

Específicamente, necesito obtener la siguiente fórmula:

y = A + Bx + e

A es la intersección y B es la pendiente, mientras que e es el error residual por punto.

Answer 1

A juzgar por el enlace que proporcionó, y mi comprensión de su problema, desea calcular la línea de mejor ajuste para un conjunto de puntos de datos. También quieres hacer esto desde los primeros principios. Esto requerirá algo de cálculo básico así como algo de álgebra lineal para resolver un sistema de ecuaciones de 2 x 2. Si recuerda de la teoría de regresión lineal, deseamos encontrar la mejor pendiente m intersección b tal que para un conjunto de puntos ([x_1,y_1], [x_2,y_2], ..., [x_n,y_n]) ( es decir, tenemos n puntos de datos), queremos minimizar la suma de residuos cuadrados entre esta línea y los puntos de datos.

En otras palabras, deseamos minimizar la función de costo F(m,b,x,y) :

b son nuestra pendiente e intersección para esta línea de mejor ajuste, mientras que y son un vector de coordenadas y que forman nuestro conjunto de datos.

Esta función es convexa, por lo que hay un mínimo óptimo que podemos determinar. El mínimo se puede determinar encontrando la derivada con respecto a cada parámetro, y estableciéndolos igual a 0. Entonces resolvemos para m y b . La intuición detrás de esto es que estamos encontrando simultáneamente myb de manera que la función de costo se minimiza conjuntamente por estos dos parámetros. En otras palabras:

OK, entonces busquemos la primera cantidad :

Podemos eliminar el factor 2 de la derivada ya que el otro lado de la ecuación es igual a 0, y también podemos hacer una distribución de términos multiplicando el término -x_i a lo largo de:

A continuación, abordemos el siguiente parámetro :

Podemos volver a soltar el factor de 2 y distribuir el -1 a lo largo de la expresión:

Sabiendo que es simplemente n , podemos simplificar lo anterior para:

Ahora, necesitamos resolver simultáneamente para b con las dos ecuaciones anteriores. Esto minimizará conjuntamente la función de costo que encuentra la mejor línea de ajuste para nuestros puntos de datos.

Al hacer una nueva reorganización, podemos aislar b en un lado de las ecuaciones y el resto en los otros lados:

Como puede ver, podemos formular esto en un sistema de ecuaciones de 2 x 2 para resolver m y b . Específicamente, reorganicemos las dos ecuaciones de arriba para que esté en forma de matriz:

Con respecto a lo anterior, podemos descomponer el problema resolviendo un sistema lineal: Ax = b . Todo lo que tienes que hacer es resolver x , que es x = A^{-1}*b . Para encontrar el inverso de un sistema 2 x 2, dada la matriz:

Lo inverso es simplemente:

Por lo tanto, al sustituir nuestras cantidades en la ecuación anterior, resolvemos para m y b en forma de matriz, y se simplifica a esto:

Llevar a cabo esta multiplicación y resolver para b individualmente, esto da:

Como tal, para encontrar la mejor pendiente e intersección que mejor se ajuste a sus datos, necesita calcular b usando las ecuaciones anteriores.

Teniendo en cuenta sus datos especificados en el enlace en sus comentarios, podemos hacer esto con bastante facilidad:

%// Define points
X = 1:5;
Y = [1 2 1.3 3.75 2.25];

%// Get total number of points
n = numel(X);

% // Define relevant quantities for finding quantities
sumxi = sum(X);
sumyi = sum(Y);
sumxiyi = sum(X.*Y);
sumxi2 = sum(X.^2);
sumyi2 = sum(Y.^2);

%// Determine slope and intercept
m = (sumxi * sumyi - n*sumxiyi) / (sumxi^2 - n*sumxi2);
b = (sumxiyi * sumxi - sumyi * sumxi2) / (sumxi^2 - n*sumxi2);

%// Display them
disp([m b])

... y obtenemos:

0.4250  0.7850

Por lo tanto, la línea de mejor ajuste que minimiza el error es:

y = 0.4250*x + 0.7850

Sin embargo, si desea usar las herramientas integradas de MATLAB, puede usar polyfit (el crédito va para Luis Mendo por proporcionar la pista). polyfit determina la línea (o n ^° orden de la curva polinómica en lugar ...) de mejor ajuste por regresión lineal al minimizar la suma de los errores al cuadrado entre la mejor línea de ajuste y sus puntos de datos. Cómo llamas a la función es así:

coeff = polyfit(x,y,order);

y son los puntos y de sus datos mientras que el order determina el orden de la línea de mejor ajuste que desea. Como ejemplo, order=1 significa que la línea es lineal, order=2 significa que la línea es cuadrática y así sucesivamente. Esencialmente, polyfit ajusta a un polinomio de orden de order dados sus puntos de datos. Dado su problema, order=1 . Como tal, dados los datos en el enlace, simplemente harías:

X = 1:5;
Y = [1 2 1.3 3.75 2.25];
coeff = polyfit(X,Y,1)

coeff =

    0.4250    0.7850

El modo en que coeff funciona es que estos son los coeficientes de la línea de regresión, comenzando desde el orden más alto en valor decreciente. Como tal, la variable coeff anterior significa que la línea de regresión se ajustó como:

y = 0.4250*x + 0.7850

El primer coeficiente es la pendiente, mientras que el segundo coeficiente es la intersección. También verá que esto coincide con el enlace que proporcionó.

Si desea una representación visual, aquí hay un gráfico de los puntos de datos, así como la línea de regresión que mejor se ajusta a estos puntos:

plot(X, Y, ''r.'', X, polyval(coeff, X));

Aquí está la trama:

polyval toma una matriz de coeficientes (normalmente producidos por polyfit ), y usted proporciona un conjunto de coordenadas x y calcula qué valores y reciben valores de x . Básicamente, estás evaluando cuáles son los puntos a lo largo de la línea de mejor ajuste.

Editar - Extendiéndose a órdenes superiores

Si desea extender para que encuentre el mejor ajuste para cualquier polinomio de orden n, no entraré en detalles, pero todo se reduce a la construcción del siguiente sistema lineal. Dada la relación para el i ^ésimo punto entre (x_i, y_i) :

Usted construiría el siguiente sistema lineal:

Básicamente, crearía un vector de puntos y , y construiría una matriz X tal que cada columna denota tomar su vector de puntos x y aplicar una operación de potencia a cada columna. Específicamente, la primera columna es la potencia cero, la primera columna es la primera potencia, la segunda columna es la segunda potencia y así sucesivamente. Lo harías hasta m , que es el polinomio de orden que deseas. El vector de e sería el error residual para cada punto en su conjunto.

Específicamente, la formulación del problema se puede escribir en forma de matriz como:

Una vez que construyes esta matriz, encontrarás los parámetros por mínimos cuadrados calculando el pseudoinverso . Cómo se deriva el pseudoinverse, puedes leerlo en el artículo de Wikipedia al que me he vinculado, pero esta es la base para minimizar un sistema por mínimos cuadrados. El pseudoinverse es la columna vertebral detrás de la minimización de mínimos cuadrados. Específicamente:

(X^{T}*X)^{-1}*X^{T} es el pseudo-inverso. X sí es una matriz muy popular, que se conoce como la matriz de Vandermonde y MATLAB tiene un comando llamado vander para ayudarte a calcular esa matriz. Una pequeña nota es que el vander en MATLAB se devuelve en orden inverso. Las potencias disminuyen de m-1 a 0. Si quieres tener esto invertido, necesitas llamar a fliplr en esa matriz de salida. Además, deberá agregar una columna más al final, que es el vector con todos sus elementos elevados a la m ^enésima potencia.

No entraré en cómo repetirías tu ejemplo para cualquier orden más alta que lineal. Te lo dejo como ejercicio de aprendizaje, pero simplemente construye el vector y , la matriz X con vander , luego encuentra los parámetros aplicando el pseudoinverso de X con el de arriba para resolver tus parámetros.

¡Buena suerte!