寻找最佳单调曲线拟合

Question

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

---

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

我最初一直在对我的结果进行 polyfit，并且在我找到一个特定的数据集之前一直运行良好。此数据集中数据的 polyfit 是非单调的。

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

使用 lsqlin。将一阶导数约束为非负 感兴趣的域的末端。

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

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

也试试这个非线性函数 F(x)。你一起使用它 lsqcurvefit 但它需要对参数进行猜测。但它是 一个很好的分析表达式，作为半经验公式给出 论文或报告。

%单调函数 F(x)，具有 c0,c1,c2,c3 变分常数 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) （看数学），但函数 % 看起来很像 x^3 %导数f(x)，概率密度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')

对于指定的输入拟合曲线：

实际上这段代码比你要求的更通用，因为多项式的次数也可以改变。

编辑：在附加点强制执行衍生约束

cmets 中指出的问题是由于派生检查仅在数据点中执行。在这些之间不执行任何检查。以下是缓解此问题的解决方案。想法：通过使用惩罚项将问题转换为无约束优化。

请注意，它使用术语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 带有非线性约束可能是更好的选择。 我让你制定细节并调整算法。我希望它足够了。