Matlabで多項式フィッティングとスプラインフィッティング

www.youtube.com

polyfit()とpolyval()を使う

clear; close all; clc;
x = [0 0.5 1.1 1.7 2.1 2.5 2.9 3.3 3.7 4.2 4.9 5.3];
y = [1.1 1.6 2.4 3.0 4.3 4.7 4.8 5.5 6.1 6.3 7.1 7.1];

for i=1:6
subplot(2,3,i);
pcoeff=polyfit(x,y,i)
xn=0:0.1:x(end);
yn=polyval(pcoeff,xn);
plot(x,y,'o');
hold on;
plot(xn,yn,'LineWidth',2);
title(sprintf('Polynomial fitting: N=%d',i));
end

N=12点に対して、n-1次（＝11次）多項式で近似した場合

clear; close all; clc;
x = [0 0.5 1.1 1.7 2.1 2.5 2.9 3.3 3.7 4.2 4.9 5.3];
y = [1.1 1.6 2.4 3.0 4.3 4.7 4.8 5.5 6.1 6.3 7.1 7.1];
n = length(x);
pcoeff=polyfit(x,y,n-1)
xn=0:0.1:x(end);
yn=polyval(pcoeff,xn);
plot(x,y,'o');
hold on;
plot(xn,yn,'LineWidth',2);

端点で値がオーバーシュートしている→Polynomial wiggleとよばれ、多項式フィッティングの欠点。

スプラインを試してみる。Matlabのinterp1()関数を使う。
www.mathworks.com
interp1()では補間メソッドを選べる。linear, nearest, next, previous, pchip, cubic, v5cubic, makima, splineがある。

clear; close all; clc;
x = [0 0.5 1.1 1.7 2.1 2.5 2.9 3.3 3.7 4.2 4.9 5.3];
y = [1.1 1.6 2.4 3.0 4.3 4.7 4.8 5.5 6.1 6.3 7.1 7.1];


subplot(231);
x2=0:0.1:x(end);
y2=interp1(x,y,x2);
plot(x,y,'o');
hold on;
plot(x2,y2,'LineWidth',2);
title('(default) linear');

subplot(232);
x2=0:0.1:x(end);
y2=interp1(x,y,x2,'nearest');
plot(x,y,'o');
hold on;
plot(x2,y2,'LineWidth',2);
title('nearest');

subplot(233);
x2=0:0.1:x(end);
y2=interp1(x,y,x2,'spline');
plot(x,y,'o');
hold on;
plot(x2,y2,'LineWidth',2);
title('spline');

subplot(234);
x2=0:0.1:x(end);
y2=interp1(x,y,x2,'previous');
plot(x,y,'o');
hold on;
plot(x2,y2,'LineWidth',2);
title('previous (=zero order hold)');

subplot(235);
x2=0:0.1:x(end);
y2=interp1(x,y,x2,'pchip');
plot(x,y,'o');
hold on;
plot(x2,y2,'LineWidth',2);
title('pchip');

subplot(236);
x2=0:0.1:x(end);
y2=interp1(x,y,x2,'makima');
plot(x,y,'o');
hold on;
plot(x2,y2,'LineWidth',2);
title('makima');

デフォルトの線形補間は普通の折れ線グラフ。Previousは0次ホールド。splineとpchipとmakimaの挙動は似ている。この3者だけ抜き出して重ねてみる。

clear; close all; clc;
x = [0 0.5 1.1 1.7 2.1 2.5 2.9 3.3 3.7 4.2 4.9 5.3];
y = [1.1 1.6 2.4 3.0 4.3 4.7 4.8 5.5 6.1 6.3 7.1 7.1];

x2=0:0.1:x(end);
y2=interp1(x,y,x2,'spline');
plot(x,y,'o');
hold on;
plot(x2,y2,'LineWidth',2);
hold on;
y3=interp1(x,y,x2,'pchip');
plot(x2,y3,'LineWidth',2);
y4=interp1(x,y,x2,'makima');
plot(x2,y4,'LineWidth',2);
legend('data','spline','pchip','makima');

Splineが微妙にはみ出してるようにもみえるけど、ほとんど同じ。