寻找最佳单调曲线拟合

Question

编辑：在我问这个问题一段时间后，一个名为MonoPoly（可在此处获取）的 R 包出现，它正是我想要的。我强烈推荐它。

---

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

我最初一直在对我的结果进行多重拟合，直到我找到一个特定的数据集之前，这一直工作得很好。该数据集中数据的多重拟合是非单调的。

我做了一些研究并在这篇文章中找到了可能的解决方案：

使用lsqlin。将感兴趣域两端的一阶导数限制为非负。

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

进一步的研究发现这篇文章推荐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')

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

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

编辑：在附加点中强制执行导数约束

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

请注意，它使用一个术语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具有非线性约束可能是更好的选择。我让你弄清楚细节并调整算法。我希望这已经足够了。