Autor Tema: Ajuste mínimos cuadrados

0 Usuarios y 1 Visitante están viendo este tema.

25 Mayo, 2020, 12:58 pm
Leído 225 veces

Bobby Fischer

  • $$\Large \color{red}\pi\,\pi\,\pi\,\pi\,\pi$$
  • Mensajes: 558
  • País: es
  • Karma: +1/-0
  • Sexo: Masculino
    • chess.com
Hola,

Aunque Matlab ya tiene la función 'polyfit', he escrito la mía propia para resolver un problema de cuadrados mínimos de tipo discreto, por si alguien lo quiere probar o aconsejarme algo.

Código: [Seleccionar]
function []=a49(x,y,n)
% El polinomio de ajuste es de grado <= n
vx=zeros(length(x),1);
vy=vx;
for i=1:length(x)
    vx(i)=x(i);
    vy(i)=y(i);
end
x=vx;
y=vy;

A=zeros(n+1);
b=zeros(n+1,1);
for i=0:n
    for j=i:n
        A(i+1,j+1)=(x.^i)'*x.^j;
    end
    b(i+1)=(x.^i)'*y;
end
A=A+triu(A,1)';
u=A\b;

ccs=zeros(n+1,1);
for j=1:n+1
    ccs(j)=u(n+2-j);
end
disp(ccs)
x1=linspace(min(x),max(x),1e3);
y1=polyval(ccs,x1);

close(figure(1))
figure(1)
hold on
plot(x,y,'rp')
plot(x1,y1,'b')
axis equal

% x=[2 4 6 8]; y=[2 11 28 40]; n=1;
% x=rand(4,1); y=rand(4,1); a49(x,y,length(x-1))


25 Mayo, 2020, 04:15 pm
Respuesta #1

mathtruco

  • Moderador Global
  • Mensajes: 4,994
  • País: cl
  • Karma: +0/-0
  • Sexo: Masculino
  • El gran profesor inspira
Hola Bobby Fischer.

No revisé el código completo, pero como pides consejos te doy el siguiente: evita al máximo los ciclos en matlab. Es bien sabido que funcionan muy lento.

Por ejemplo,

Código: [Seleccionar]
vx=zeros(length(x),1);
vy=vx;

for i=1:length(x)
    vx(i)=x(i);
    vy(i)=y(i);
end

usa

Código: [Seleccionar]
vx=[x ; 0]; % teniendo cuidado si es traspuesta o no
vy=[y ; 0];

Creo que el otro ciclo también puedes evitarlo.

25 Mayo, 2020, 07:51 pm
Respuesta #2

Abdulai

  • Moderador Global
  • Mensajes: 2,387
  • País: ar
  • Karma: +0/-0
  • Sexo: Masculino
Hola,
Aunque Matlab ya tiene la función 'polyfit', he escrito la mía propia para resolver un problema de cuadrados mínimos de tipo discreto, por si alguien lo quiere probar o aconsejarme algo.
Si buscás opinión uno opina.

Es el armado de las matrices que puede abreviarse y acelerarse, pues como comenta mathtruco, matlab es lento en los bucles.

Código: [Seleccionar]
vx=zeros(length(x),1);
vy=vx;
for i=1:length(x)
    vx(i)=x(i);
    vy(i)=y(i);
end
x=vx;
y=vy;

puede escribirse como
Código: [Seleccionar]
if size(x,1)==1 ; x=x' ; end
if size(y,1)==1 ; y=y' ; end

y
Código: [Seleccionar]
A=zeros(n+1);
b=zeros(n+1,1);
for i=0:n
    for j=i:n
        A(i+1,j+1)=(x.^i)'*x.^j;
    end
    b(i+1)=(x.^i)'*y;
end
A=A+triu(A,1)';

por
Código: [Seleccionar]
A = x.^[0:n] ;
b= A'*y ;
A = A'*A ;

o bien, la solución del sistema como
Código: [Seleccionar]
A = x.^[0:n] ;
u = pinv(A)*y ;

25 Mayo, 2020, 09:12 pm
Respuesta #3

Bobby Fischer

  • $$\Large \color{red}\pi\,\pi\,\pi\,\pi\,\pi$$
  • Mensajes: 558
  • País: es
  • Karma: +1/-0
  • Sexo: Masculino
    • chess.com
Hola,

Gracias por vuestra ayuda, lo tendré en cuenta.

Saludos.

25 Mayo, 2020, 09:18 pm
Respuesta #4

Abdulai

  • Moderador Global
  • Mensajes: 2,387
  • País: ar
  • Karma: +0/-0
  • Sexo: Masculino
Ah!   La inversión del vector tampoco necesita bucle.

Código: [Seleccionar]
ccs=zeros(n+1,1);
for j=1:n+1
    ccs(j)=u(n+2-j);
end

es  lo mismo que
Código: [Seleccionar]
ccs = flip(u) ;

Sacando los bucles y las variables auxiliares la resolución del sistema te queda
Código: [Seleccionar]
if size(x,1)==1 ; x=x' ; end
if size(y,1)==1 ; y=y' ; end

ccs = flip(pinv(x.^[0:n])*y) ;

O sin flip :)    ccs = pinv(x.^[n:-1:0])*y ;

25 Mayo, 2020, 11:37 pm
Respuesta #5

Bobby Fischer

  • $$\Large \color{red}\pi\,\pi\,\pi\,\pi\,\pi$$
  • Mensajes: 558
  • País: es
  • Karma: +1/-0
  • Sexo: Masculino
    • chess.com
Siempre tengo en mente que lo que hago se puede hacer mejor, pero no es hasta que alguien te lo muestra que efectivamente ves que el listón puede subirse.

Gracias.