# Aproximación por mínimos cuadrados

## Aproximación por mínimos cuadrados continua

Aproximamos la función 1-x^3 en [0,1] mediante un polinomio de grado 1 o menor,
con el producto escalar usual entre funciones (el dado por la integral del producto de las funciones). 

f(x):=1-x^3;

Intentamos crear un array con las funciones que forman una base del espacio de polinomios considerado.

g(x):=[1,x];

g(3.0);

g(x)[1];

g(x)[2];

basesize:length(g(x));

Introducimos la matriz del sistema y el término independiente.

msistema:zeromatrix(basesize,basesize);

terindep:zeromatrix(basesize,1);

for i :1 thru basesize  do(
  for j:1 thru basesize do(

  msistema[i,j]: integrate(g(x)[i]*g(x)[j], x, 0, 1)
),
terindep[i,1]: integrate(g(x)[i]*f(x), x, 0, 1)
);

La matriz de coeficientes también se puede declarar usando la función auxiliar

h[i,j]:= integrate(g(x)[j]*g(x)[i],x,0,1)

y utilizando un genmatrix(). De forma análoga para el término independiente.

print(msistema,".y=",terindep);

Ahora ya tenemos definida la matriz del sistema y el término independiente. 
Lo resolvemos con alguno de nuestros métodos anteriores.

A : copymatrix(msistema) $
b : copymatrix(terindep) $
n : basesize $
tol : 10^(-15)$
for k:1 while k < n do (
    filabuena : 0,
    for i:k while i <= n do (
        if A[i,k] > tol then (
            filabuena : i,
            i : n+1
        )
    ),
    if filabuena = 0 then error("ERROR: Matriz singular"),
    if filabuena > k then (
        print("cero en ", k, ", pivote en ", filabuena),
        tmp : A[k],
        A[k] : A[filabuena],
        A[filabuena] : tmp,
        tmp : b[k],
        b[k] : b[filabuena],
        b[filabuena] : tmp
    ),
    for i:k+1 while i<=n do (
        m : -A[i,k]/A[k,k],
        A[i] : A[i] + m*A[k],
        b[i] : b[i] + m*b[k]
    )
) $
print(A,".y=",b) $

 y: zeromatrix (n, 1)$
for i:n thru 1 step -1 do 
y[i] : (b[i] - sum(A[i,j]*y[j],j,i+1,n))/A[i,i] $
print("y=",y) $

aprox(x):=sum(x^(l-1)*y[l,1],l,1,basesize);

aprox(x);

wxplot2d([f(x),aprox(x)],[x,0,1]);

Calculamos el error cuadrático medio:

integrate((f(x)- aprox(x))^2,x,0,1);

## Aproximación por mínimos cuadrados discreta

Aproximamos los puntos {(1,2),(2,3),(3,5),(4,7),(5,3)} mediante un polinomio de grado dos o menor.

abcisas:[1,2,3,4,5];

ordenadas:[2,3,5,7,3];

deg:2;

Declaramos la matriz auxiliar y con ella la matriz del sistema y el término independiente.

eme:zeromatrix(length(abcisas),deg+1);

for i:1 thru length(abcisas) do(
  for j:1 thru deg+1 do(
    eme[i,j]:abcisas[i]^(j-1)
)
);

eme;

msist:transpose(eme).eme;

termindep:transpose(eme).transpose(ordenadas);

dim:deg+1;

Resolvemos el sitema lineal asociado

A : copymatrix(msist) $
b : copymatrix(termindep) $
n : dim $
tol : 10^(-15)$
for k:1 while k < n do (
    filabuena : 0,
    for i:k while i <= n do (
        if A[i,k] > tol then (
            filabuena : i,
            i : n+1
        )
    ),
    if filabuena = 0 then error("ERROR: Matriz singular"),
    if filabuena > k then (
        print("cero en ", k, ", pivote en ", filabuena),
        tmp : A[k],
        A[k] : A[filabuena],
        A[filabuena] : tmp,
        tmp : b[k],
        b[k] : b[filabuena],
        b[filabuena] : tmp
    ),
    for i:k+1 while i<=n do (
        m : -A[i,k]/A[k,k],
        A[i] : A[i] + m*A[k],
        b[i] : b[i] + m*b[k]
    )
) $
print(A,".y=",b) $

y: zeromatrix (n, 1)$
for i:n thru 1 step -1 do 
y[i] : (b[i] - sum(A[i,j]*y[j],j,i+1,n))/A[i,i] $
print("y=",y) $

Representamos el ajuste

xy : create_list([abcisas[i],ordenadas[i]],i,1,length(abcisas));

apro:sum(y[l,1]*x^(l-1),l,1,dim);

wxplot2d( [[discrete,xy], apro], [x, 0, 6], 
        [style, [points,3,2], [lines,1,1]]) $