最良の単調曲線フィットを見つける

Question

編集：私はこの質問をした後で、MonoPoly（利用可能なhere）と呼ばれるRパッケージが出てきました。私は非常にそれをお勧めします。

---

カーブに合わせたい点があります。曲線は単調でなければならず（値が決して減少しない）、すなわち、曲線は上に行くか、または平らにしかならない。

私はもともと私の結果をポリフィッティングしていましたが、これは特定のデータセットが見つかるまで素晴らしい結果を出していました。このデータセットのデータのポリフィットは単調ではありませんでした。

私はいくつかの研究を行なったし、thisポストで可能な解決策が見つかりました：

使用lsqlinを。関心のあるドメインの両端の の両方で、一次導関数を非負に制限する。

私は数学のバックグラウンドではなく、プログラミングから来ているので、これは私を少し超えています。彼が言ったように、最初の派生物を非負であるように制約する方法はわかりません。また、私は私の場合、私はlsqcurvefitを使用する必要があるので、私は曲線が必要だと思うが、私は単調曲線を生成するためにそれを制約する方法を知らない。

さらなる研究がlsqcurvefitを推薦thisポストを上げたが、私は重要な部分を使用する方法を見つけ出すことはできません。

も、この非線形関数F（X）を試してみてください。 lsqcurvefitと一緒に使用しますが、パラメータの推測が必要です。しかし、それは の論文や報告書で半経験的な式として与えられる素晴らしい分析式です。

%の単調関数F（x）は、C0、C1、C2、C3 varitional定数を有するF（X）= C3 + EXP（C0 - C1^2 /（4 * C2））SQRT（PI） 【数1】Erfi（x）= erf（i * x）/（2 * sqrt（c2）））/（2 * sqrt（c2））x（x）= dF/dx = exp（c0 + c1 * x + 3）％（x + 3）％の確率密度関数f（x） c2 * x。^ 2）

私は単調なカーブを持つ必要がありますが、どのようにするのかはわかりません。この情報。乱数が「推測開始」に十分であるかどうか。 lsqcurvefitは最高ですか？どのようにして最良のフィット感のある単調曲線を生成するには？ここで

おかげ

Answer 1

はlsqlinを使った簡単なソリューションです。微分制約は各データポイントに適用され、必要に応じて簡単に変更できます。

2つの係数行列が必要です。最小二乗誤差計算では1つ（C）、データポイントでは派生要素では1つ（A）です。指定された入力の場合

% Following lsqlin's notations 

%-------------------------------------------------------------------------- 
% PRE-PROCESSING 
%-------------------------------------------------------------------------- 
% for reproducibility 
rng(125) 
degree = 3; 
n_data = 10; 
% dummy data 
x  = rand(n_data,1); 
d  = rand(n_data,1) + linspace(0,1,n_data).'; 

% limit on derivative - in each data point 
b  = zeros(n_data,1); 

% coefficient matrix 
C  = nan(n_data, degree+1); 
% derivative coefficient matrix 
A  = nan(n_data, degree); 

% loop over polynomial terms 
for ii = 1:degree+1 
    C(:,ii) = x.^(ii-1); 
    A(:,ii) = (ii-1)*x.^(ii-2); 
end 

%-------------------------------------------------------------------------- 
% FIT - LSQ 
%-------------------------------------------------------------------------- 
% Unconstrained 
% p1 = pinv(C)*y 
p1 = fliplr((C\d).') 

p2 = polyfit(x,d,degree) 

% Constrained 
p3 = fliplr(lsqlin(C,d,-A,b).') 

%-------------------------------------------------------------------------- 
% PLOT 
%-------------------------------------------------------------------------- 
xx = linspace(0,1,100); 

plot(x, d, 'x') 
hold on 
plot(xx, polyval(p1, xx)) 
plot(xx, polyval(p2, xx),'--') 
plot(xx, polyval(p3, xx)) 

legend('data', 'lsq-pseudo-inv', 'lsq-polyfit', 'lsq-constrained', 'Location', 'southoutside') 
xlabel('X') 
ylabel('Y') 

近似曲線：

多項式の次数は、同様に変更することができますので、実際にこのコードは、あなたが要求したものよりもより一般的です。

編集：問題誘導体チェックがデータポイントにのみ適用されることに起因しているコメントで指摘見所

で誘導体制約を強制します。その間にチェックは行われません。以下は、この問題を緩和するための解決策です。アイデア：ペナルティ項を使用して問題を制約のない最適化に変換します。

penという用語を使用して、デリバティブチェックの違反をペナルティするため、真の最小二乗誤差解決法ではないことに注意してください。さらに、結果はペナルティ関数に依存します。

function lsqfit_constr 
% Following lsqlin's notations 

%-------------------------------------------------------------------------- 
% PRE-PROCESSING 
%-------------------------------------------------------------------------- 
% for reproducibility 
rng(125) 
degree = 3; 

% data from comment 
x  = [0.2096 -3.5761 -0.6252 -3.7951 -3.3525 -3.7001 -3.7086 -3.5907].'; 
d  = [95.7750 94.9917 90.8417 62.6917 95.4250 89.2417 89.4333 82.0250].'; 
n_data = length(d); 

% number of equally spaced points to enforce the derivative 
n_deriv = 20; 
xd  = linspace(min(x), max(x), n_deriv); 

% limit on derivative - in each data point 
b  = zeros(n_deriv,1); 

% coefficient matrix 
C  = nan(n_data, degree+1); 
% derivative coefficient matrix 
A  = nan(n_deriv, degree); 

% loop over polynom terms 
for ii = 1:degree+1 
    C(:,ii) = x.^(ii-1); 
    A(:,ii) = (ii-1)*xd.^(ii-2); 
end 

%-------------------------------------------------------------------------- 
% FIT - LSQ 
%-------------------------------------------------------------------------- 
% Unconstrained 
% p1 = pinv(C)*y 
p1  = (C\d); 
lsqe = sum((C*p1 - d).^2); 

p2  = polyfit(x,d,degree); 

% Constrained 
[p3, fval] = fminunc(@error_fun, p1); 

% correct format for polyval 
p1  = fliplr(p1.') 
p2 
p3  = fliplr(p3.') 
fval 

%-------------------------------------------------------------------------- 
% PLOT 
%-------------------------------------------------------------------------- 
xx  = linspace(-4,1,100); 

plot(x, d, 'x') 
hold on 
plot(xx, polyval(p1, xx)) 
plot(xx, polyval(p2, xx),'--') 
plot(xx, polyval(p3, xx)) 

% legend('data', 'lsq-pseudo-inv', 'lsq-polyfit', 'lsq-constrained', 'Location', 'southoutside') 
xlabel('X') 
ylabel('Y') 

%-------------------------------------------------------------------------- 
% NESTED FUNCTION 
%-------------------------------------------------------------------------- 
    function e = error_fun(p) 
     % squared error 
     sqe = sum((C*p - d).^2); 
     der = A*p; 

     % penalty term - it is crucial to fine tune it 
     pen = -sum(der(der<0))*10*lsqe; 

     e = sqe + pen; 
    end 
end 

グラデーション自由方法は、例えば、正確に派生制約を強制することによって、問題を解決するために使用される可能性があります：非線形制約付き

[p3, fval] = fminsearch(@error_fun, p_ini); 

%-------------------------------------------------------------------------- 
% NESTED FUNCTION 
%-------------------------------------------------------------------------- 
function e = error_fun(p) 
    % squared error 
    sqe = sum((C*p - d).^2); 
    der = A*p; 

    if any(der<0) 
     pen = Inf; 
    else 
     pen = 0; 
    end 

    e = sqe + pen; 
end 

fminconより良い選択かもしれません。 詳細を解説し、アルゴリズムを調整することができます。私はそれが十分であることを願っています。