MATLAB Examples

function fit_test(which)
% This function tests the fit_nl and fit_nl_ex routines. Several fits in 1D-3D
% with built-in functions or custom defined functions are performed.
%
% Syntax
%   function fit_test()
%   function fit_test(which)
%
% Input parameter
%   which   which test 1-7 (or nothing for a menu)

if nargin == 0

    % option which not given, print a menu
    which = 1;
    while which ~= 0
        % the menu
        fprintf('Which test would You like to do?\n');
        fprintf(' (1) Calling built-in function in fit_func()\n');
        fprintf(' (2) 1D fit of noisy data with a Lorentzian function\n');
        fprintf(' (3) 2D fit of noisy data with a custom defined function\n');
        fprintf(' (4) 3D fit of noisy data with a Gaussian function\n');
        fprintf(' (5) 1D fit of noisy data with 3 Lorentzian peaks\n');
        fprintf(' (6) 2D fit of noisy data with 4 Gaussian peaks and complex side conditions\n');
        fprintf(' (7) Runtime comparison of lsqcurvefit(), fmincon() and fminsearch()\n');
        fprintf(' (Enter) Exit\n');
        which = input('Choose [Enter]:', 's');
        if isempty(which)
            % which is empty when enter was pressed
            which = 0;
        else
            which = int32(str2double(which));
            if which > 0 && which < 8
                % call the right test
                fit_test(which);
            else
                %  all other cases (not only Enter) we also exit
                which = 0;
            end
        end
    end
else

    % according to which execute a test
    fprintf('\n');
    switch which
        case 1
            % a test: Calling built-in functions in fit_func() ------------

            disp('Calling built-in function in fit_func()');

            % x values run from - 1 to 5
            xdata.x = -1:0.1:5;

            % function parameter for function 1, a gaussian with amplitude
            % 10, center at 0,5, FWHM 2 and background 0.1
            x = [0.1, 10, 0.5, 2];
            % compute it using a built-in function
            y1 = fit_func('1d_gaussian', x, xdata);

            % using the same x-values we have the parameters for function
            % 2, a sine with amplitude 4, cycle duration 3, start at -1
            % and background of 3
            x = [3., 4., 3, -1];
            y2 = fit_func(@custom_func_1dsine, x, xdata);

            % plot both in one graph
            figure;
            hold on;
            plot(xdata.x, y1, 'b');
            plot(xdata.x, y2, 'r');
            title('function data created by a call to fit\_func()');

        case 2
            % tests 2-4 use function fit_nl()
            % a test: 1D fitting with a Lorentzian ------------------------

            disp('1D fit of noisy data with a Lorentzian function');

            % prepare image and initial parameters
            xdata.x = 0.5:99.5; % 100 x-values between 0 and 100

            % the 1D Lorentzian model has 4 parameters: background,
            % amplitude, center and FWHM
            x_true = [0.5, 20., 35., 10.];

            % left and right borders
            lb   = [0., 1., 0., 0.1];
            ub   = [1., 100., 100., 20.];

            % calculate the function corresponding to x_true
            y_true = fit_func('1d_lorentzian', x_true, xdata);

            % add gaussian noise to y_true to obtain y_meas
            y_meas = random('Normal', y_true, 3); % standard deviation = 3

            % set reasonable start values for fitting parameters
            x0 = [0.1, 10., 50., 5.];

            % fit y_meas with start values x0 with a call to fit_nl
            % parameter fixed (5th) is not set
            x_fit = fit_nl('1d_lorentzian', x0, xdata, y_meas, [], lb, ub);

            % calculate y_fit corresponding to x_fit
            y_fit = fit_func('1d_lorentzian', x_fit, xdata);

            % plot everything in a figure
            figure;
            hold on;
            plot(xdata.x, y_true, 'b');
            plot(xdata.x, y_meas, 'k');
            plot(xdata.x, y_fit, 'r');
            ylim([-10 25]);
            title('test 1 - 1D fitting of a Lorentzian peak');
            text(2, -8, 'blue - original curve, black - noisy curve, red - fitted curve');

            % display parameters overview
            disp('original parameters');
            disp(x_true);
            disp('fitted parameters');
            disp(x_fit);

        case 3
            % test 2D fitting of a sine with a custom curve ---------------

            % Comment: Here we use a custom function and give two fields x
            % and y in xdata so that the function knows on which grid in
            % the R^2 (the 2D world) the function should be computed. x and
            % y have the same size as the output and can be conveniently by
            % created by a call to ndgrid(). However, every other way,
            % where xdata is a string, a cell array or whatsever, which
            % makes your function know whose x and y values (of the
            % independent variables) to take is also okay.

            disp('2D fit of noisy data with a custom defined function');

            % prepare xdata / computational grid
            [x, y]  = ndgrid(0.5:19.5, 0.5:19.5);
            xdata.x = x;
            xdata.y = y;
            % our custom function has 4 parameters (see custom_func_2dsine
            % for a function expression)

            % the true parameter
            x_true = [10., 0.1, 0.1, 0.2];

            % calculate the 2D function corresponding to x_true
            y_true = fit_func(@custom_func_2dsine, x_true, xdata);

            % add gaussian noise to y_true to obtain y_meas
            y_meas = random('Normal', y_true, 6); % standard deviation = 6

            % set reasonable start values for fitting parameters
            x0 = [8., 0.12, 0.12, 0];

            % fit y_meas with start values x0 with a call to fit_nl
            % parameter fixed, lb and ub are not set
            x_fit = fit_nl(@custom_func_2dsine, x0, xdata, y_meas);

            % calculate y_fit corresponding to x_fit
            y_fit = fit_func(@custom_func_2dsine, x_fit, xdata);

            % plot everything in a figure
            clims = [-100 220];
            subplot2(2, 2, 1);
            imagesc(y_true, clims);
            title('original image');
            subplot2(2, 2, 2);
            imagesc(y_meas, clims);
            title('noisy image');
            subplot2(2, 2, 3);
            imagesc(y_fit, clims);
            title('fitted image');
            subplot2(2, 2, 4);
            imagesc(y_true - y_fit, clims);
            title('difference original - fitted image');

            % display parameters overview
            disp('original parameters');
            disp(x_true);
            disp('fitted parameters');
            disp(x_fit);

        case 4
            % test 3D fitting of a gaussian -------------------------------

            disp('3D fit of noisy data with a Gaussian function');

            % prepare xdata / computational grid
            [x, y, z]  = ndgrid(0.5:19.5, 0.5:19.5, 0.5:19.5);
            % in all directions 20 points ranging vom 0.5 to 19.5
            xdata.x = x;
            xdata.y = y;
            xdata.z = z;

            % the 3D Gaussian function has 7 parameters (see fit_func.m)

            % the true parameters
            x_true = [0.1, 50., 8., 11., 9., 2., 1.5, 3.];

            % calculate the 2D function corresponding to x_true
            y_true = fit_func('3d_gaussian', x_true, xdata);

            % add gaussian noise to y_true to obtain y_meas
            y_meas = random('Normal', y_true, 6); % standard deviation = 6

            % set reasonable start values for fitting parameters
            x0 = [0.5, 47., 7., 12., 8.5, 3., 1.8, 2.7];

            % the first parameter (background) shall be kept fixed
            fixed = [1, 0, 0, 0, 0, 0, 0, 0];

            % fit y_meas with start values x0 with a call to fit_nl
            % parameter lb and ub are not set
            [x_fit resnorm] = fit_nl('3d_gaussian', x0, xdata, y_meas, fixed);

            % calculate y_fit corresponding to x_fit
            y_fit = fit_func('3d_gaussian', x_fit, xdata);

            % we do not plot (would need to be a 3D plot)
            fprintf('3D fit completed with fixed parameter 1, residual norm is %.2f\n', resnorm);

            % display parameters overview
            disp('original parameters');
            disp(x_true);
            disp('initial guess for parameters');
            disp(x0);
            disp('fitted parameters');
            disp(x_fit);

        case 5
            % test 5 and 6 use function fit_nl_ex()
            % test: 1D fit of 3 Lorentzian peaks --------------------------

            disp('1D fit of noisy data with 3 Lorentzian peaks');

            % prepare xdata / computational grid
            xdata = [];
            xdata.x  = 0.5:99.5;

            % Three 1D Lorentzian functions have 1+3*(2*1+1)=10 parameter
            % (see fitfunc.m)

            % the true parameters
            x_true = [0.1, 20., 20., 8., 10., 50., 5., 15., 70., 6.5];

            % calculate the 2D function corresponding to x_true
            y_true = fit_func('1d_lorentzian', x_true, xdata);

            % add gaussian noise to y_true to obtain y_meas
            y_meas = random('Normal', y_true, 3); % standard deviation = 3

            % set reasonable start values for fitting parameters
            x0 = [0.3, 18., 23., 7., 12., 55., 3., 17., 75., 4.5];

            % create options
            options.opt = optimset('Display', 'iter');

            % fit y_meas with start values x0 with a call to fit_nl_ex
            x_fit = fit_nl_ex('1d_lorentzian', x0, xdata, y_meas, options);

            % since we did not do any special, the fit did exactly the same
            % as with fit_nl, but this demonstrates only, that fmincon and
            % lsqcurvefit relies on the same principles

            % calculate y_fit corresponding to x_fit
            y_fit = fit_func('1d_lorentzian', x_fit, xdata);

            % plot everything in a figure
            figure;
            hold on;
            plot(xdata.x, y_true, 'b');
            plot(xdata.x, y_meas, 'k');
            plot(xdata.x, y_fit, 'r');
            ylim([-10 25]);
            title('test 5 - 1D fitting of 3 Lorentzian peaks');
            text(2, -4, 'blue - original curve, black - noisy curve, red - fitted curve');
            text(2, -6, 'watch the center peak which is sometimes not fitted correctly');
            text(2, -8, 'because of noise and a bad starting point');

            % display parameters overview
            disp('original parameters');
            disp(x_true);
            disp('fitted parameters');
            disp(x_fit);

        case 6
            % test: 2D of 4 gaussian peaks --------------------------------

            disp('2D fit of noisy data with 4 Gaussian peaks and complex side conditions');
            % now our showcase :)))

            % prepare xdata / computational grid
            xdata = [];
            [x, y] = ndgrid(0.5:99.5, 0.5:99.5);
            xdata.x = x;
            xdata.y = y;

            % 4 2D Gaussian functions have 1+4*(2*2+1)=21 parameter
            % thats a lot (see fitfunc.m)

            % the true parameters
            x_true = 0.1;   % background
            x_true = [x_true, 20., 30., 30., 8., 10.];  % first peak
            x_true = [x_true, 22., 70., 60., 6., 9.];   % second peak
            x_true = [x_true, 26., 40., 70., 11., 7.];  % third peak
            x_true = [x_true, 20., 80., 90., 8., 8.];   % fourd peak

            % calculate the 2D function corresponding to x_true
            y_true = fit_func('2d_gaussian', x_true, xdata);

            % add gaussian noise to y_true to obtain y_meas
            y_meas = random('Poisson', y_true); % standard deviation = sqrt(mean)

            % set reasonable start values for fitting parameters
            x0 = 0.3; % background
            x0 = [x0, 15., 25., 35., 8.5, 8.5];     % first peak
            x0 = [x0, 20., 65., 65., 8.5, 8.5];     % second peak
            x0 = [x0, 30., 45., 65., 8.5, 8.5];     % third peak
            x0 = [x0, 15., 75., 85., 8.5, 8.5];     % fourd peak

            % calculate the 2D function corresponding to x0
            y0 = fit_func('2d_gaussian', x0, xdata);

            % create options
            options.lb = zeros(size(x0)); % all parameters should be at least positive
            options.ub = zeros(size(x0)) + 100; % and at most 100
            options.likelihood  = 'Poisson';
            options.penalty_fun = @penalty_fun;

            % fit y_meas with start values x0 with a call to fit_nl_ex
            fprintf('Fitting in progress, please wait!');
            x_fit = fit_nl_ex('2d_gaussian', x0, xdata, y_meas, options);

            % compute image corresponding to x_fit
            y_fit = fit_func('2d_gaussian', x_fit, xdata);

            % plot everything in a figure
            clims = [0 30];
            subplot2(2, 2, 1);
            imagesc(y_true, clims);
            colormap(hot);
            title('original image');
            subplot2(2, 2, 2);
            imagesc(y_meas, clims);
            title('noisy image');
            subplot2(2, 2, 3);
            imagesc(y0, clims);
            title('initial guess image');
            subplot2(2, 2, 4);
            imagesc(y_fit, clims);
            title('fitted image');

            % display parameters overview
            fprintf('parameters\n');
            fprintf('true\t initial\t fitted\n');
            for ki = 1:length(x_true)
                fprintf('%.1f\t %.1f\t %.1f\n', x_true(ki), x0(ki), x_fit(ki));
            end

        case 7
            % test: Runtime comparison lsqcurvefit/fmincon/fminsearch -----

            disp('Measure running time for a simple example: 1D Gaussian peak');
            % prepare a simple example - a 1D gaussian peak

            % All default options for maxIter, TolX, etc. should be the same
            % for lsqcurvefit(), fmincon(), fminsearch()
            % A good sign for this is that the reached residual values were
            % the same in my case

            % prepare image and initial parameters
            xdata.x = 0.5:99.5; % 100 x-values between 0 and 100

            % the 1D Gaussian model has 4 parameters: background,
            % amplitude, center and FWHM
            x_true = [0.5, 20., 35., 10.];

            % calculate the function corresponding to x_true
            y_true = fit_func('1d_gaussian', x_true, xdata);

            % to have more runs to add we create several noisy copies of
            % our data and fit them all
            N = 100;
            fprintf('Performing %d fitting starts per functions.\n', N);

            % add gaussian noise to y_true to obtain y_meas
            y_true2 = repmat(y_true, N, 1);
            y_meas2 = random('Normal', y_true2, 3); % standard deviation = 3

            % set reasonable start values for fitting parameters
            x0 = [0.1, 15., 40., 8.];

            % call to lsqcurvefit via fit_nl
            fprintf('Fitting via lsqcurvefit\n');
            tic
            r1 = 0;
            for ki = 1 : N
                [x r] = fit_nl('1d_gaussian', x0, xdata, y_meas2(ki, :));
                r1 = r1 + r;
            end
            r1 = r1 / N;
            toc
            fprintf('Mean residual value (resnorm of lsqcurvefit) %.2f\n', r1);

            % call to fmincon via fit_nl_ex
            fprintf('Fitting via fmincon\n');
            tic
            r2 = 0;
            for ki = 1 : N
                [x r] = fit_nl_ex('1d_gaussian', x0, xdata, y_meas2(ki, :));
                r2 = r2 + r;
            end
            r2 = r2 / N;
            toc
            fprintf('Mean residual value (fval of fmincon) %.2f\n', r2);

            % call to fminsearch via fit_nl_ex
            options = [];
            options.use_fmincon = 'no';
            fprintf('Fitting via fminsearch\n');
            tic
            r3 = 0;
            for ki = 1 : N
                [x r] = fit_nl_ex('1d_gaussian', x0, xdata, y_meas2(ki, :), options);
                r3 = r3 + r;
            end
            r3 = r3 / N;
            toc
            fprintf('Mean residual value (fval of fminsearch) %.2f\n', r3);

        otherwise
            error('Variable which out of range!');
    end
    fprintf('\n');

end

% local functions ---------------------------------------------------------

% custom function for test 1
    function y = custom_func_1dsine(x, xdata)
        % has 4 parameters and expects xdata.x to hold the x values
        % form: p1 + p2 * sin(((x-p4)/p3)*2pi)
        y = zeros(size(xdata.x)) + x(1);
        y = y + x(2) * sin((xdata.x - x(4)) / x(3) * 2 * pi);
    end

% custom function for test 3
    function y = custom_func_2dsine(x, xdata)
        % is 2D (expects xdata.x and xdata.y)
        % form: p1 * sin(p2 * x) * cos(p3 * y) + p4
        y = x(1) * sin(x(2) * xdata.x) * cos(x(3) * xdata.y) + x(4);
    end

% penalty function for test 6
    function y = penalty_fun(x)
        % simply adds a quadratic penalty if amplitude of the 4 gaussian peaks
        % is far away from 22, e.g. this would help against overfitting with
        % darker (relative to true brightness) peaks if the number of peaks
        % would be unknown, a weighting factor should be introduced
        % additionally
        y = (x(2)-22)^2+(x(7)-22)^2+(x(12)-22)^2+(x(17)-22)^2;
    end

end
Which test would You like to do?
 (1) Calling built-in function in fit_func()
 (2) 1D fit of noisy data with a Lorentzian function
 (3) 2D fit of noisy data with a custom defined function
 (4) 3D fit of noisy data with a Gaussian function
 (5) 1D fit of noisy data with 3 Lorentzian peaks
 (6) 2D fit of noisy data with 4 Gaussian peaks and complex side conditions
 (7) Runtime comparison of lsqcurvefit(), fmincon() and fminsearch()
 (Enter) Exit
Error using ==> input
Cannot call INPUT from EVALC.

Error in ==> fit_test at 27
        which = input('Choose [Enter]:', 's');