[Ricardo Boza]

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))



[mathtruco]

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.

[Abdulai]

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 ;


[Ricardo Boza]

Hola,

Gracias por vuestra ayuda, lo tendré en cuenta.

Saludos.

[Abdulai]

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 ;

[Ricardo Boza]

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.