matlab - 寻找最佳单调曲线拟合

Question

编辑:在我问这个问题一段时间后，一个名为 MonoPoly (可用 here )的 R 包出现了，它完全符合我的要求。我强烈推荐它。

我有一组要拟合曲线的点。曲线必须是单调的(值永远不会减少)，即曲线只能向上或保持平坦。

我最初一直在对我的结果进行多拟合，并且在我找到一个特定的数据集之前一直很好用。该数据集中数据的 polyfit 是非单调的。

我做了一些研究，并在 this 帖子中找到了一个可能的解决方案:

Use lsqlin. Constrain the first derivative to be non-negative at both ends of the domain of interest.

我来自编程而不是数学背景，所以这有点超出我的范围。正如他所说，我不知道如何将一阶导数限制为非负。另外，我认为在我的情况下我需要一条曲线，所以我应该使用 lsqcurvefit 但我不知道如何约束它以产生单调曲线。

进一步的研究出现了 this 推荐 lsqcurvefit 的帖子，但我不知道如何使用重要的部分:

Try this non-linear function F(x) also. You use it together with lsqcurvefit but it require a start guess on the parameters. But it is a nice analytic expression to give as a semi-empirical formula in a paper or a report.

%Monotone function F(x), with c0,c1,c2,c3 varitional constants F(x)= c3 + exp(c0 - c1^2/(4*c2))sqrt(pi)... Erfi((c1 + 2*c2*x)/(2*sqrt(c2))))/(2*sqrt(c2))

%Erfi(x)=erf(i*x) (look mathematica) but the function %looks much like x^3 %derivative f(x), probability density f(x)>=0 f(x)=dF/dx=exp(c0+c1*x+c2*x.^2)

我必须有一个单调曲线，但我不知道如何去做，即使有所有这些信息。随机数是否足以进行“开始猜测”。 lsqcurvefit 是最好的吗？我如何使用它来生成最佳拟合单调曲线？

谢谢

Answer 1

这是一个使用 lsqlin 的简单解决方案。在每个数据点强制执行导数约束，如果需要，可以轻松修改。

需要两个系数矩阵，一个 ( C ) 用于最小二乘误差计算，另一个 ( 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 可能是更好的选择。
我让你解决细节并调整算法。我希望它足够了。