使用单个句柄存储和访问分段插值

Question

我有一个包含两个点向量的数据集，X 和 Y 表示 "exponential-like" 现象（即 Y = A*exp(b*x)）的测量值。当用指数方程拟合它时，我得到了一个 漂亮的 拟合，但是当我用它来计算东西时，事实证明拟合不如我准确希望。

目前我最有前途的想法是分段指数拟合（每次取大约 6 (x,y) 对），在我手动测试的情况下它似乎提供了更好的结果。这是一个图表来说明我的意思：

// Assuming a window of size WS=4:
 - - - - - - - - - - - -  //the entire X span of the data
[- - - -]- - // the fit that should be evaluated for X(1)<= x < X(floor(WS/2)+1)
 -[- - - -]- // the fit that should be evaluated for X(3)<= x < X(4)
...
 - - - - - -[- - - -]- - // X(8)<= x < X(9)
... //and so on

一些注意事项：

我考虑过在拟合之前过滤数据，但这很棘手，因为我对它包含的噪声类型一无所知。
我希望使用单个函数句柄可以访问分段拟合（包括所有不同的情况）。它与 MATLAB 的 Shape-preserving "PCHIP" interpolant 非常相似，只是它应该使用指数函数。
拟合的创建不需要在另一个代码的运行期间发生。我什至更喜欢单独创建它。
我不担心最终函数的潜在不平滑。
我能想到的解决这个问题的唯一方法是定义一个 .m 文件，如 Fit a Curve Defined by a File 中所述，但这将需要手动编写条件，几乎与我要点一样多（或者明明写个代码帮我生成这段代码，也是相当费力的）。

相关代码片段：

clear variables; clc;
%% // Definitions
CONSTS.N_PARAMETERS_IN_MODEL = 2; %// For the model y = A*exp(B*x);
CONSTS.WINDOW_SIZE = 4; 
assert(~mod(CONSTS.WINDOW_SIZE,2) && CONSTS.WINDOW_SIZE>0,...
    'WINDOW_SIZE should be a natural even number.');
%% // Example Data
X = [0.002723682418630,0.002687088539567,0.002634005004610,0.002582978173834,...
     0.002530684550171,0.002462144527884,0.002397219225698,0.002341097974950,...
     0.002287544321171,0.002238889510803]';
Y = [0.178923076435990,0.213320004768074,0.263918364216839,0.324208349386613,...
     0.394340431220509,0.511466688684182,0.671285738221314,0.843849959919278,...
     1.070756756433334,1.292800046096531]';
assert(CONSTS.WINDOW_SIZE <= length(X),...
    'WINDOW_SIZE cannot be larger than the amount of points.');
X = flipud(X); Y = flipud(Y); % ascending sorting is needed later for histc
%% // Initialization
nFits = length(X)-CONSTS.WINDOW_SIZE+1;
coeffMat(nFits,CONSTS.N_PARAMETERS_IN_MODEL) = 0; %// Preallocation
ft = fittype( 'exp1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
%% // Populate coefficient matrix
for ind1 = 1:nFits
    [xData, yData] = prepareCurveData(...
        X(ind1:ind1+CONSTS.WINDOW_SIZE-1),Y(ind1:ind1+CONSTS.WINDOW_SIZE-1));
    %// Fit model to data:
    fitresult = fit( xData, yData, ft, opts );
    %// Save coefficients:
    coeffMat(ind1,:) = coeffvalues(fitresult);
end
clear ft opts ind1 xData yData fitresult ans
%% // Get the transition points
xTrans = [-inf; X(CONSTS.WINDOW_SIZE/2+1:end-CONSTS.WINDOW_SIZE/2); inf];

此时，xTrans 和 coeffMat 包含评估拟合所需的所有信息。为了证明这一点，我们可以查看一些测试数据的向量：

testPts = [X(1), X(1)/2, mean(X(4:5)), X(CONSTS.WINDOW_SIZE)*1.01,2*X(end)];

...最后应该在句柄内部大致发生以下情况：

%% // Get the correct fit# to be evaluated:
if ~isempty(xTrans) %// The case of nFits==1
    rightFit = find(histc(testPts(3),xTrans));
else
    rightFit = 1;
end
%% // Actually evaluate the right fit
f = @(x,A,B)A*exp(B*x);
yy = f(testPts(3),coeffMat(rightFit,1),coeffMat(rightFit,2));

所以我的问题是如何在接受任意大小的单个句柄内保存最后一位（连同所有拟合系数）输入 个要插值的点？

相关资源：

whosebug.com/questions/16777921/matlab-curve-fitting-exponential-vs-linear/

Answer 1

还不是很清楚，但为什么不把东西放到 class 中呢？

classdef Piecewise < handle

    methods

        % Construction
        function [this] = Piecewise(xmeas, ymeas)
            ... here compute xTrans and coeffMat...
        end

        % Interpolation
        function [yinterp] = Evaluate(xinterp)
            ... Here use previously computed xTrans and coeffMat ...
        end
    end

    properties(SetAccess=Private, GetAccess=Private)
        xTrans;
        coeffMat;
    end

end

通过这种方式，您可以在构造期间预先计算 xTrans 向量和 coeffMat 矩阵，然后在需要评估 [=18= 中的 xinterp 值处的插值时重用这些属性]方法。

% The real measured data
xmeas = ...
ymeas = ...

% Build piecewise interpolation object
piecewise = Piecewise(x,meas, ymeas);

% Rebuild curve at any new positions
xinterp = ...
yinterp = piecewise.Evaluate(xinterp);

函数式语法

如果你真的喜欢函数句柄这样的语法，你仍然可以在内部使用上面的Piecewise对象并添加额外的static 方法到 return 它作为函数句柄：

classdef Piecewise < handle

    ... see code above...

    methods(Static=true)
        function [f] = MakeHandle(xmeas, ymeas)
        %[
            obj = Piecewise(xmeas, ymeas);
            f = @(x)obj.Evaluate(x);
        %]
    end

end

可以这样使用：

f = Piecewise.MakeHandle(xmeas, ymeas);
yinterp = f(xinterp);

PS1：如果您绝对想强制使用此语法，则稍后可以将 Evaluate 和 Piecewise 构造函数方法设为私有。

PS2：要完全隐藏面向对象的设计，您可以将 MakeHandle 变成完全 classic 例程（将与静态和用户不必在 MakeHandle).

前面输入 Piecewise.

没有oop的最后一个解决方案

将所有内容放在一个 .m 文件中：

function [f] = MakeHandle(xmeas, ymeas)
    ... Here compute xTrans and coeffMat ...
    f = @(x)compute(x, xTrans, coeffMat);% Passing xTrans/coeffMatt as hidden parameters
end
function [yinterp] = compute(xinterp, xTrans, coeffMat)
    ... do interpolation here...
 end

Answer 2

作为 CitizenInsane 答案的扩展，以下是允许“handle-y”访问内部 Evaluate 函数的替代方法。

    function b = subsref(this,s)
       switch s(1).type
          case '()'
             xval = s.subs{:};
             b = this.Evaluate(xval);
          otherwise %// Fallback to the default behavior for '.' and '{}':
             builtin('subsref',this,s);
       end
    end

参考文献：docs1, docs2, docs3, question1

P.S。 docs2 被引用是因为我最初尝试做 subsref@handle （这是调用超类方法，正如在 OOP 中覆盖子类中的方法时所期望的那样），但这在 MATLAB 中不起作用，它相反需要 builtin() 来实现相同的功能。

使用单个句柄存储和访问分段插值

Storing and accessing a piecewise interpolant using a single handle

matlab

interpolation

handle

curve-fitting

piecewise

一些注意事项：

相关代码片段：

相关资源：