非线性最小二乘法曲线拟合

具体细节参见代码：

function [beta] = nonLinearSquareFit(X,Y,beta,option)
% X,Y 要拟合的数值对，一维向量/矩阵
% beta - 参数初始值以及输出的参数值，一维向量/矩阵
% option - 优化参数

if nargin < 4
    option.tolFun = 1.0e-8;  % 阈值
    option.tolX = 1.0e-8;    % 阈值
    option.DerivStep = eps^(1/3);  % 
    option.MaxIter = 200;  % 最大迭代次数
end
% X与Y的元素数量必须一致
if numel(X)~=numel(Y)
    error('the length of X must be equal to Y!')
end
X = reshape(X,[numel(X), 1]);
y = reshape(Y,[numel(Y), 1]);
beta = reshape(beta, [numel(beta), 1]);

betatol = option.tolX;
rtol = option.tolFun;
DerivStep = option.DerivStep;
if isscalar(DerivStep)
    DerivStep = repmat(DerivStep, size(beta));
else
    DerivStep = reshape(DerivStep, size(beta));
end
maxiter = option.MaxIter;

% 更新参数的系数
lambda = 0.01; 
sqrteps = sqrt(eps(class(beta)));

yfit = model(X, beta);
% 判断yfit和y是否有NAN
nans = ((isnan(y(:))) | isnan(yfit(:)));
r = y(:) - yfit(:);
r(nans) = [];
p = numel(beta);
sse = r'*r;

zerosp = zeros(p,1,class(r));
iter = 0;

while iter < maxiter
    iter = iter + 1;
    betaold = beta;
    sseold = sse;
    
    J = getjcb(X, betaold, DerivStep, yfit, nans);
    diagJtJ = sum(abs(J).^2, 1);
    Jplus = [J; diag(sqrt(lambda*diagJtJ))];
    rplus = [r; zerosp];
%     Jplus = [J];
%     rplus = [r];
    step = Jplus \ rplus;
%     step1 = solvingEquat(Jplus, rplus);
%     step2 = pinv(Jplus) * rplus;
%     fprintf(2, 'diff step1 = %f\n', step(:) - step1(:));
%     fprintf(2, 'diff step2 = %f\n', step(:) - step2(:));
    beta(:) = beta(:) + step(:);
    
    yfit = model(X, beta);
    fullr = y(:) - yfit(:);
    r = fullr(~nans);
    sse = r'*r;
    
    if sse < sseold
        lambda = max(0.1*lambda,eps);
    else
        while sse > sseold
            lambda = 10*lambda;
            if lambda > 1e16
                break
            end
            Jplus = [J; diag(sqrt(lambda*sum(J.^2,1)))];
            step = Jplus \ rplus;
            beta(:) = betaold(:) + step;
            yfit = model(X, beta);
            fullr = y(:) - yfit(:);
            r = fullr(~nans);
            sse = r'*r;
        end
    end
    
    if norm(step) < betatol*(sqrteps+norm(beta))
        break
    elseif abs(sse-sseold) <= rtol*sse
        break
    end
end

function [J] = getjcb(X, beta, DerivStep, yfit, nans)
    % get jocabian
    beta = reshape(beta, [1, numel(beta)]);
    [numThetaRows, numParams] = size(beta);
    
    delta = zeros(numThetaRows, numParams, 'double');
    J = zeros(numel(yfit), numParams, 'double');
    for i = 1:numParams
        delta(:,i) = DerivStep(i) * beta(:,i);
        deltaZero = (delta(:,i) == 0);
        if any(deltaZero)
            nTheta = sqrt(sum(beta(deltaZero,:).^2, 2));
            delta(deltaZero,i) = DerivStep(i) * (nTheta + (nTheta==0));
        end
        thetaNew = beta + delta;
        yplus = model(X, thetaNew);
        yplus(nans) = [];
        dy = yplus(:) - yfit(:);
        J(:,i) = dy./delta(1,i);
        delta(:,i) = 0;
    end
    


function [Y] = model(X,beta)
    % 待拟合的曲线
    Y = beta(1)*exp(-X/beta(2)) + beta(3);

测试代码如下：

%% name: testNonLinearSquare.m
%% test nlinfit
    X = [100,144,167,200,300,400,600,1000];
    Y = [103,100,98,97,94,93,91,90];

    func = inline('beta(1)*exp(-x/beta(2))+beta(3)','beta','x');
    beta0 = [10 144 1]';  % 初始参数值，其值关系到拟合的效果
    beta1 = nlinfit(X,Y,func,beta0);

    [beta2] = nonLinearSquareFit(X,Y,beta0);

    xi = (50 : 1100);
    yi1 = beta1(1)*exp(-xi/beta1(2)) +beta1(3);
    yi2 = beta2(1)*exp(-xi/beta2(2)) +beta2(3);
    diff = sum(yi2-yi1)
    
    figure
    plot(X,Y,'.', xi,yi1,'r-', xi,yi2,'b-');
    radStr{1} = 'y vs x data1';
    radStr{2} = 'nlinfit';
    radStr{3} = 'nonLinearSquareFit';
    legend(radStr);
    xlabel('X');
    ylabel('Y');
    title('curve fitting');

输出的拟合曲线为：