Aproximación de funciones por mínimos cuadrados

Ejemplo: Una empresa ha estimado los costos unitarios de un producto (en miles de pesos) en función de la partida de producción medida en miles de unidades:

Producción
(miles)
Costo unitario
(miles de $)
1.42.5
1.53
24
2.55
3.55.5

A partir de las siguientes propuestas de aproximación se pide determinar, con la que mejor represente los datos, el costo de producción de 3000 unidades

\[ f_{1}(x) = c_{1}x^{2} + c_{2}x + c_{3} \]

\[ f_{2}(x) = c_{1}x^{2} + c_{2}e^{x} \]


Resolución por mínimos cuadrados

Función \(f_1(x) = c_1 x^2 + c_2 x + c_3\)

Paso 1 (\(f_1\)): Planteo del sistema de ecuaciones normales

Funciones base: \(\phi_1(x)=x^2\), \(\phi_2(x)=x\), \(\phi_3(x)=1\). El sistema es:

\[ \begin{bmatrix} \sum x^2\cdot x^2 & \sum x\cdot x^2 & \sum 1\cdot x^2 \\ \sum x^2\cdot x & \sum x\cdot x & \sum 1\cdot x \\ \sum x^2\cdot 1 & \sum x\cdot 1 & \sum 1\cdot 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \end{bmatrix} = \begin{bmatrix} \sum y\,x^2 \\ \sum y\,x \\ \sum y \end{bmatrix} \]

Tabla de cálculos auxiliares:

\(x\)\(y\)\(x^2\)\(x^3\)\(x^4\)\(x^2 y\)\(x y\)\(y\)
1.42.51.962.7443.84164.93.52.5
1.53.02.253.3755.06256.754.53.0
2.04.04.008.00016.00016.08.04.0
2.55.06.2515.62539.062531.2512.55.0
3.55.512.2542.875150.062567.37519.255.5
\(\sum\)2026.7172.619214.0291126.27547.7520

El sistema resulta:

\[ \begin{cases} 214.0291\,c_1 + 72.619\,c_2 + 26.71\,c_3 = 126.275 \\ 72.619\,c_1 + 26.71\,c_2 + 10.9\,c_3 = 47.75 \\ 26.71\,c_1 + 10.9\,c_2 + 5\,c_3 = 20 \end{cases} \]
Paso 2 (\(f_1\)): Resolución por Gauss (con pivoteo parcial)

El método de Gauss consiste en transformar el sistema en una matriz triangular superior (ceros debajo de la diagonal) y luego resolver de abajo hacia arriba (sustitución hacia atrás). El pivoteo parcial agrega un paso extra: antes de eliminar cada columna, se busca en esa columna (de la fila actual para abajo) el elemento con mayor valor absoluto y se lo lleva a la posición de pivote intercambiando filas completas. Esto se hace para evitar dividir por números muy chicos, lo cual amplifica errores de redondeo.

Notación: en la matriz aumentada, cada celda se identifica como \(a_{\text{fila},\text{columna}}\) y el término independiente (lo que está a la derecha de la barra) como \(b_{\text{fila}}\). Cada fila completa se llama \(F_1, F_2, F_3\):

\[ \left[\begin{array}{ccc|c} a_{11} & a_{12} & a_{13} & b_1 \\ a_{21} & a_{22} & a_{23} & b_2 \\ a_{31} & a_{32} & a_{33} & b_3 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \\ & \leftarrow F_3 \end{aligned} \]

Matriz aumentada inicial (el sistema del paso 1 puesto en forma de matriz):

\[ \left[\begin{array}{ccc|c} 214.0291 & 72.619 & 26.71 & 126.275 \\ 72.619 & 26.71 & 10.9 & 47.75 \\ 26.71 & 10.9 & 5 & 20 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \\ & \leftarrow F_3 \end{aligned} \]

Eliminación columna 1 — Pivoteo parcial:

Se quiere hacer cero todo lo que está debajo de \(a_{11}\) en la columna 1. Pero primero, pivoteo: se busca cuál de los tres valores de la columna 1 tiene el mayor valor absoluto, y si no está en \(F_1\), se intercambian las filas completas para ponerlo ahí.

  • \(|a_{11}| = 214.0291\) ← es el mayor, ya está en \(F_1\)
  • \(|a_{21}| = 72.619\)
  • \(|a_{31}| = 26.71\)

Como 214.0291 ya está en la fila 1 (que es donde necesitamos el pivote), no hay intercambio. Si, por ejemplo, el mayor hubiera sido 72.619 (en \(F_2\)), habría que intercambiar \(F_1 \leftrightarrow F_2\) completas antes de continuar.

Ahora se calculan los multiplicadores: son los factores por los que hay que multiplicar \(F_1\) para restársela a cada fila de abajo y así hacer cero en la columna 1. Se obtienen dividiendo cada elemento de la columna 1 por el pivote \(a_{11}\):

\[ m_{21} = \frac{a_{21}}{a_{11}} = \frac{72.619}{214.0291} = 0.3393 \]

\[ F_2 = F_2 - m_{21} \cdot F_1 \] \[ a'_{22} = 26.71 - 0.3393 \times 72.619 = 2.0707 \] \[ a'_{23} = 10.9 - 0.3393 \times 26.71 = 1.8374 \] \[ b'_2 = 47.75 - 0.3393 \times 126.275 = 4.9055 \] \[ \left[\begin{array}{ccc|c} 214.0291 & 72.619 & 26.71 & 126.275 \\ 0 & 2.0707 & 1.8374 & 4.9055 \\ 26.71 & 10.9 & 5 & 20 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \\ & \leftarrow F_3 \end{aligned} \]


\[ m_{31} = \frac{a_{31}}{a_{11}} = \frac{26.71}{214.0291} = 0.1248 \] \[ F_3 = F_3 - m_{31} \cdot F_1 \] \[ a'_{32} = 10.9 - 0.1248 \times 72.619 = 1.8374 \] \[ a'_{33} = 5 - 0.1248 \times 26.71 = 1.6667 \] \[ b'_3 = 20 - 0.1248 \times 126.275 = 4.2414 \] \[ \left[\begin{array}{ccc|c} 214.0291 & 72.619 & 26.71 & 126.275 \\ 0 & 2.0707 & 1.8374 & 4.9055 \\ 0 & 1.8374 & 1.6667 & 4.2414 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \\ & \leftarrow F_3 \end{aligned} \]
\[ m_{32} = \frac{a'_{32}}{a'_{22}} = \frac{1.8374}{2.0707} = 0.8873 \] \[ F_3 = F_3 - m_{32} \cdot F_2 \] \[ a''_{33} = 1.6667 - 0.8873 \times 1.8374 = 0.0363 \] \[ b''_3 = 4.2414 - 0.8873 \times 4.9055 = -0.1115 \]

Matriz triangular superior (objetivo cumplido: todos ceros debajo de la diagonal):

\[ \left[\begin{array}{ccc|c} 214.0291 & 72.619 & 26.71 & 126.275 \\ 0 & 2.0707 & 1.8374 & 4.9055 \\ 0 & 0 & 0.0363 & -0.1115 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \\ & \leftarrow F_3 \end{aligned} \]

Sustitución hacia atrás:

Se resuelve de abajo hacia arriba. La última fila tiene una sola incógnita, la penúltima dos, etc.

De \(F_3\):   \(0.0363 \cdot c_3 = -0.1115\)

\[ c_3 = \frac{-0.1115}{0.0363} = -3.0717 \]

De \(F_2\):   \(2.0707 \cdot c_2 + 1.8374 \cdot c_3 = 4.9055\). Ya se conoce \(c_3\), se reemplaza y se despeja \(c_2\):

\[ c_2 = \frac{4.9055 - 1.8374 \times (-3.0717)}{2.0707} = \frac{4.9055 + 5.6435}{2.0707} = 5.0946 \]

De \(F_1\):   \(214.0291 \cdot c_1 + 72.619 \cdot c_2 + 26.71 \cdot c_3 = 126.275\). Ya se conocen \(c_2\) y \(c_3\):

\[ c_1 = \frac{126.275 - 72.619 \times 5.0946 - 26.71 \times (-3.0717)}{214.0291} = \frac{126.275 - 370.047 + 82.024}{214.0291} = -0.7552 \]

Resultado:

\[ \color{limegreen} \boxed{ \color{black} f_1(x) = -0.7552\,x^2 + 5.0946\,x - 3.0717 } \]

Función \(f_2(x) = c_1 x^2 + c_2 e^x\)

Paso 1 (\(f_2\)): Planteo del sistema de ecuaciones normales

Funciones base: \(\phi_1(x)=x^2\), \(\phi_2(x)=e^x\). El sistema es:

\[ \begin{bmatrix} \sum x^2\cdot x^2 & \sum e^x\cdot x^2 \\ \sum x^2\cdot e^x & \sum e^x\cdot e^x \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} \sum y\,x^2 \\ \sum y\,e^x \end{bmatrix} \]

Tabla de cálculos auxiliares (con \(e^x\) aproximado):

\(x\)\(y\)\(x^2\)\(e^x\)\(x^4\)\(x^2 e^x\)\(e^{2x}\)\(x^2 y\)\(e^x y\)
1.42.51.964.05523.84167.94816.4454.910.138
1.53.02.254.48175.062510.08420.0856.7513.445
2.04.04.007.389116.00029.55654.59816.029.556
2.55.06.2512.182539.062576.141148.4131.2560.912
3.55.512.2533.1155150.0625405.661096.667.375182.14
\(\sum\)2026.7161.224214.0291529.3891336.138126.275296.191

Sistema resultante:

\[ \begin{bmatrix} 214.0291 & 529.389 \\ 529.389 & 1336.138 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} 126.275 \\ 296.191 \end{bmatrix} \]
Paso 2 (\(f_2\)): Resolución por Gauss (con pivoteo parcial)

Notación: misma convención que en \(f_1\). Ahora es una matriz 2×2 (columna 1 = coeficiente de \(c_1\), columna 2 = coeficiente de \(c_2\)):

\[ \left[\begin{array}{cc|c} a_{11} & a_{12} & b_1 \\ a_{21} & a_{22} & b_2 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \end{aligned} \]

Matriz aumentada inicial:

\[ \left[\begin{array}{cc|c} 214.0291 & 529.393 & 126.275 \\ 529.393 & 1336.175 & 296.187 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \end{aligned} \]

Eliminación columna 1 — Pivoteo parcial:

Se busca el mayor valor absoluto en la columna 1:

  • \(|a_{11}| = 214.0291\)
  • \(|a_{21}| = 529.393\) ← es el mayor, está en \(F_2\)

Como el mayor no está en \(F_1\), se intercambian las filas completas \(F_1 \leftrightarrow F_2\). Notar que el 1336.175 viaja junto con su fila (no se eligió, solo acompaña al 529.393):

\[ \left[\begin{array}{cc|c} 529.393 & 1336.175 & 296.187 \\ 214.0291 & 529.393 & 126.275 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \end{aligned} \]

Multiplicador:

\[ m_{21} = \frac{a_{21}}{a_{11}} = \frac{214.0291}{529.393} = 0.4043 \]

\(F_2 \leftarrow F_2 - m_{21} \cdot F_1\) (a cada elemento de \(F_2\) se le resta \(m_{21}\) veces el elemento correspondiente de \(F_1\)):

\[ a'_{22} = 529.393 - 0.4043 \times 1336.175 = -10.811 \] \[ b'_2 = 126.275 - 0.4043 \times 296.187 = 6.529 \]

Matriz triangular superior:

\[ \left[\begin{array}{ccc|c} 529.393 & 1336.175 & c_2 & 296.187 \\ 0 & -10.811 & c_1 & 6.529 \end{array}\right] \quad \begin{aligned} & \leftarrow F_1 \\ & \leftarrow F_2 \end{aligned} \]

Sustitución hacia atrás:

De \(F_2\):   \(-10.811 \cdot c_2 = 6.529\)

\[ c_2 = \frac{6.529}{-10.811} = -0.6039 \]

De \(F_1\):   \(529.393 \cdot c_1 + 1336.175 \cdot c_2 = 296.187\). Ya se conoce \(c_2\):

\[ c_1 = \frac{296.187 - 1336.175 \times (-0.6039)}{529.393} = \frac{296.187 + 807.100}{529.393} = 2.0838 \]

Resultado:

\[ \color{limegreen} \boxed{ \color{black} f_2(x) = 2.0838\,x^2 - 0.6039\,e^x } \]

Nota: la profesora resolvió con pivoteo total (intercambiando filas y columnas para poner el mayor de toda la matriz, 1336.175, en la posición (1,1)). Ambos métodos dan el mismo resultado; la diferencia es que en pivoteo total hay que recordar que al intercambiar columnas las incógnitas se reordenan.

Paso 3: Cálculo del error cuadrático de ambas funciones

Se necesita medir qué tan bien cada función se ajusta a los datos originales. Para eso se calcula el error cuadrático: para cada punto \((x_i, y_i)\) de la tabla original, se evalúa la función y se calcula la diferencia con el valor real \(y_i\). Esa diferencia se eleva al cuadrado (para que no se cancelen errores positivos con negativos) y se suman todas. Cuanto menor sea \(S\), mejor se ajusta la función a los datos.

\(x\) \(y\) \(f_1(x)\) \((y - f_1)^2\) \(f_2(x)\) \((y - f_2)^2\)
1.42.52.58050.00651.63530.7477
1.53.02.87100.01661.98211.0362
2.04.04.09670.00943.87290.0161
2.55.04.94480.00305.66670.4445
3.55.55.50820.00015.52810.0008
\(S_1\) 0.0356 \(S_2\) 2.2454
Paso 4: Comparación de errores y evaluación para 3000 unidades

Se tienen:

Claramente \(S_1 \ll S_2\), por lo tanto la función \(f_1\) es la que mejor representa los datos.

Evaluamos \(f_1\) en \(x = 3\) (3000 unidades, recordar que la producción está en miles):

\[ f_1(3) = -0.7552 \times 9 + 5.0946 \times 3 - 3.0717 \] \[ f_1(3) = -6.7968 + 15.2838 - 3.0717 \] \[ f_1(3) = 5.4153 \ (\text{miles de pesos}) \]

Por lo tanto, el costo estimado es de $5415.3.

Nota: Los coeficientes pueden variar ligeramente según la precisión usada en las sumatorias y en la resolución del sistema; el procedimiento es el fundamental.