38.2: Encontrar la mejor solución en un sistema sobredeterminado

%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import sympy as sym
import time
sym.init_printing(use_unicode=True)

Dejar\(Ax=y\) ser un sistema de ecuaciones\(m\) lineales en\(n\) variables. Una solución de mínimos cuadrados de\(Ax=y\) es una solución\(\hat{x}\) en\(R^n\) tal que:

\[ \min_{\hat{x}}\|y - A\hat{x}\| \nonumber \]

Tenga en cuenta que\(y\) sustituimos nuestra variable típica\(b\) aquí porque usaremos\(b\) más adelante para representar la intercepción a una línea y queremos intentar evitar confusiones en la notación. También concuerda con la imagen de arriba.

En otras palabras,\(\hat{x}\) es un valor de\(x\) para el cual\(Ax\) es lo más cercano posible a\(y\). De conferencias anteriores, sabemos que esto es cierto si el vector\(y-A\hat{x}\) es ortogonal (perpendicular) al espacio de columna de\(A\).

También sabemos que el producto punto es cero si dos vectores son ortogonales. Así que tenemos\(a \cdot (Ax - y) = 0 \), para todos los vectores\(a\) en el espacio de columna de\(A\).

Las columnas de\(A\) abarcan el espacio de columna de\(A\). Denotar las columnas de\(A\) as\(A = [a_1, \cdots, a_n]\). Entonces tenemos\(a_1^\top (Ax - y) = 0\), es lo mismo que tomar la transposición de\(A\) y hacer una matriz multiplicar\(A^\top (Ax - y) = 0\).

\[a_1^\top (Ax - y) = 0, \\ a_2^\top(Ax-y)=0\\\vdots \\a_n^\top(Ax-y)=0 \nonumber \]

Es decir:

\[A^\top Ax = A^\top y \nonumber \]

La ecuación anterior se llama la solución de mínimos cuadrados a la ecuación original\(Ax=y\). La matriz\(A^\top A\) es simétrica e invertible. Entonces resolver para se\(\hat{x}\) puede calcular de la siguiente manera:

\[x = (A^\top A)^{-1}A^\top y \nonumber \]

La matriz también\((A^\top A)^{-1}A^\top\) se llama inversa izquierda.

Ejemplo

Un investigador ha realizado experimentos de una dosis de Hormona particular en una población de ratas. En la tabla se muestra el número de muertes en cada nivel de dosificación probado. Determinar la línea de mínimos cuadrados y utilizarla para predecir el número de muertes de ratas a dosis de hormona de 22.

Nivel hormonal	20	25	30	35	40	45	50
Fallecidos	101	115	92	64	60	50	49

H = [20,25,30,35,40,45,50]
f = [101,115, 92,64,60,50,49]

plt.scatter(H,f)
plt.xlabel('Hormone Level')
plt.ylabel('Fatalities')
f = np.matrix(f).T

Queremos determinar una línea que se exprese mediante la siguiente ecuación

\[f = aH + b, \nonumber \]

para aproximar la conexión entre la dosis de Hormona (\(H\)) y las muertes\(f\). Es decir, queremos encontrar\(a\) (pendiente) e\(b\) (intercepción y) para esta línea. Primero definimos la variable\ (
x =\ left [
\ begin {matrix}
a\\
b
\ end {matrix}
\ right]
\) como el vector de columna que necesita ser resuelto.

Hacer esto

Reescribe el sistema de ecuaciones a la forma\(Ax=y\) definiendo tus matrices numpy A e y usando los datos de arriba:

#put your code here

Pregunta

Calcula la matriz cuadrada\(C = A^\top A\) y el vector del lado derecho modificado como\(A^\top y\) (Llámenlo Ati):

#put your code here

Pregunta

Encuentra la solución de mínimos cuadrados resolviendo\(Cx=A^\top y\) para\(x\).

# Put your code here

Pregunta

Dada la solución anterior, definir los dos escalares pendiente a e y-interceptar b.

#put your code here

El siguiente código trazará los datos originales y la línea estimada por los coeficientes encontrados en la cuación anterior.

H = [20,25,30,35,40,45,50]
f = [101,115, 92,64,60,50,49]
plt.scatter(H,f)

H2 = np.linspace(np.min(H), np.max(H))

f2 = a * H2 + b

plt.plot(H2, f2)

Pregunta

Repita el análisis anterior pero ahora con un polinomio de ocho órdenes.

Pregunta

Juega con la función interactiva a continuación ajustando el grado de aproximación de ajuste de mínimos cuadrados. Después extiende los parámetros x_min y x_max. ¿Crees que un polinomio de ocho órdenes es un buen modelo para este conjunto de datos? ¿Por qué o por qué no?

from ipywidgets import interact, fixed
import ipywidgets as widgets

@interact(x=fixed(H), y=fixed(f), degree=widgets.IntSlider(min=1, max=8, step=1, value=8), x_min=widgets.IntSlider(min=min(H)-10, max=min(H), step=1, value=min(H)), x_max=widgets.IntSlider(min=max(H), max=max(H)+10, step=1, value=max(H)))
def graphPolyN(x, y, x_min, x_max, degree):
    p = np.polyfit(x, y, degree)
    f = np.poly1d(p)
    
    x_pred = np.linspace(x_min, x_max, 1000)
    y_pred = f(x_pred)
    
    plt.scatter(x, y, color="red")
    plt.plot(x_pred, y_pred)

Pregunta

Consulta el rango de\(C=A^\top A\) para el caso anterior. ¿Qué obtienes? ¿Por qué?