フィッティングされた曲線が特定の点を通るような制約を行う場合、線形制約を与えることが必要です。 そのため、POLYFIT 関数、及び Curve Fitting Toolbox の関数では、線形制約を規定することができません。
なお、 Optimization Toolbox の LSQLIN 関数を使用することで、このような制約を与えることができます。
例えば、次のコマンドで生成されるデータについて考えます。
c = [1 -2 1 -1];
x = linspace(-2,4);
y = c(1)*x.^3+c(2)*x.^2+c(3)*x+c(4) + randn(1,100);
plot(x,y,'.b-')
次のように実行すると、 POLYFIT を使って 3 次の多項式として制約のないフィッティングの結果のプロットを表示することができます。
hold on
c = polyfit(x,y,3);
yhat = c(1)*x.^3+c(2)*x.^2+c(3)*x+c(4);
plot(x,yhat,'r','linewidth',2)
しかしながら特定の点、例えば (x0, y0) を通るよう制約を与えたい場合には、
x0 = 1;
y0 = 10;
次の記述例のように Optimization Toolbox の LSQLIN 関数を使い、線形最小二乗問題を解くことで、特定の点を通るようなフィッティングの結果を得ることができます。
x = x(:); %列ベクトルとする
y = y(:);
% "C" は "x" の Vandermonde 行列
n = 3; % 多項式の次数
V(:,n+1) = ones(length(x),1,class(x));
for j = n:-1:1
V(:,j) = x.*V(:,j+1);
end
C = V;
% "d" は目標値 "y" のベクトル
d = y;
% 今回の場合 不等式制約はない , すなわち、以下
A = [];
b = [];
% 曲線が特定の点にヒットすることを強いるために、線形制約を使う
% 今回は、 "Aeq" は "x0" の Vandermoonde 行列
Aeq = x0.^(n:-1:0);
% "beq" が、その点で曲線が取るべき値
beq = y0;
p = lsqlin( C, d, A, b, Aeq, beq )
% フィッティング結果の曲線の評価のため、 POLYVAL 関数を使う
yhat = polyval( p, x );
% オリジナルデータの表示 :
plot(x,y,'.b-')
hold on
% 特定の点の表示 :
plot(x0,y0,'gx','linewidth',4)
% フィッテイングの結果の表示 :
plot(x,yhat,'r','linewidth',2)
hold off
