fit_nl_ex

function [x, resnorm] = fit_nl_ex(model, x0, xdata, ydata, options)
% Uses fmincon or fminsearch to perform MLE parameter estimation with the
% additonal possibility to penalize the likelihood.
%
% This function uses the power of Matlab functions fmincon or fminsearch to
% minimize a constructed negative log. likelihood function which calculates
% for a given noise model (Gaussian/Poisson) the negative log. likelihood
% of parameters x0 and a model function to give a measurement ydata. It is
% therefore equivalent to fit parameters x0 of a function to data, but
% unlike lsqcurvefit() it is not restricted to least square problems.
% Additionally the possibility to use built-in functions (see fitfunc.m)
% and additional penalties which depend on the parameters in the
% likelihood. With this a-priori information about distributions of parameters
% can be incorporated.
%
% Syntax
%   function [x, resnorm] = fit_nl_ex(model, x0, xdata, ydata, options)
%
% Input parameters
%   model   Defines the function which is used for the fit. This is either
%           a keyword for built-in functions (see fit_func.m) or a function
%           pointer to a function equivalent to fun used by fmincon/fminsearch,
%           i.d. function F = myfun(x, xdata)
%   x0      The initial guess for the model parameter to be estimated.
%           For the nubmer of meaning of parameter for the built-in
%           functions, see fit_func.m
%   xdata   Input data for the model function/indepedent variables.
%           For the xdata which is required by the built-in functions see
%           fit_func.m
%   ydata   Output data to be matched by the model function
%
%   Optional parameter
%   options A struct with optional fields
%       fixed   Array of 0 or 1 determining which parameters should stay fixed
%               in the fitting process.
%               default: set to zeros(size(x0) - no fixed parameter
%               Note: Even if some parameter are set fixed the model function
%               still should take all parameter (including the fixed).
%       lb      Lower bound for x
%               default: set to -Infinity
%       ub      Upper bound for x
%               default: set to +Infinity
%       opt     the options which can be forwarded to fmincon/fminsearch
%               default: all warnings and printing will be suppressed
%       use_fmincon     ={'yes'(default),'no'} determines which function
%                       from Matlab to use
%       penalty_fun     penalty function pointer for a function of type
%                       L = fun(x) which computes an additional factor for
%                       the neg.log.likelihood depending on the actual x
%                       default: not set
%       likelihood      ={'Gaussian'(default), 'Poisson'} determines which
%                       probability noise model is used to compute the
%                       likelihood
%       A, b, Aeq, beq, nonlcon
%           Like in the call to fmincon, are ignored if fminsearch is used
%           instead.
%
% Output parameters
%   x       Resulting paramters that minimize the distance between model
%           and ydata
%   resnorm Quadratic norm of (model(x, xdata) - ydata)
%
% Comment
%   When using fmincon this function can do everything that fit_nl() can do
%   and extends this even a little bit.

% parameter checking and setting to default values ------------------------

% need at least 4 parameter
if nargin < 4
    error('Not enough arguments given!');
end

x0    = double(x0); % fmincon/fminsearch need double input
ydata = double(ydata);

s = warning('query', 'all'); % save state of warnings

if nargin < 5
    options = [];
end

% developer-comment: we cannot use
% ~isfield(options, 'fixed') | isempty(options.fixed) because Matlab will
% always evaluate also the second condition even if the first already gave
% true and therefore we need more code which we pack into a function
if isnotfieldorisempty(options, 'fixed')
    options.fixed = zeros(size(x0));
end
if isnotfieldorisempty(options, 'lb')
    options.lb = zeros(size(x0)) - Inf;
end
if isnotfieldorisempty(options, 'ub')
    options.ub = zeros(size(x0)) + Inf;
end
if isnotfieldorisempty(options, 'opt')
    options.opt = optimset('Display','off');
    warning('off', 'all');
end
if isnotfieldorisempty(options, 'use_fmincon')
    options.use_fmincon = 'yes';
end
if isnotfieldorisempty(options, 'penalty_fun')
    options.penalty_fun = 0;
    % a check with isa(.., 'function_handle') will then result in false
end
if isnotfieldorisempty(options, 'likelihood')
    options.likelihood = 'Gaussian';
end
% for options.A/b/Aeq/beq/nonlcon the default is []
if ~isfield(options, 'A')
    options.A = [];
end
if ~isfield(options, 'b')
    options.b = [];
end
if ~isfield(options, 'Aeq')
    options.Aeq = [];
end
if ~isfield(options, 'beq')
    options.beq = [];
end
if ~isfield(options, 'nonlcon')
    options.nonlcon = [];
end

if size(options.fixed) ~= size(x0) | size(options.lb) ~= size(x0) | size(options.ub) ~= size(x0)
    error('Parameter x0 and options.fixed/lb/ub have not same size!');
end

% we do not check for integers, although Poissonian noise only gives
% integers, but sometimes some devices deliver significant decimal places
% and it does not hurt the algorithm, negative values however will surely
% result in havoc
if strcmp(options.likelihood, 'Poisson')
    if any(ydata < 0)
        error('Likelihood type is Poisson but image data is negative!');
    end
end

% assemble 'global' struct accessible from subfunction F

p = [];
p.fix = options.fixed;              % fixed parameters
p.x   = x0;                         % initial parameters (needed for knowing which are fixed)
p.fun = model;                      % the model function (either keyword or function pointer)
p.lik = options.likelihood;         % the likelihood type
p.ydata = ydata;                    % the measured data
p.xdata = xdata;                    % the independent variable
p.pen_fun = options.penalty_fun;    % penalty function pointer

% reduce x0 to parameters not fixed (also lb and ub)
x0 = x0(options.fixed == 0);
options.lb = options.lb(options.fixed == 0);
options.ub = options.ub(options.fixed == 0);

% the call for the internal Matlab function
switch options.use_fmincon
    case 'yes'
        % all options are given
        [x, resnorm] = fmincon(@L, x0, options.A, options.b, options.Aeq, options.beq, options.lb, options.ub, options.nonlcon, options.opt);
    case 'no'
        % with fminsearch, some options must be ignored
        [x, resnorm] = fminsearch(@L, x0, options.opt);
    otherwise
        error('Unknown content of variable options.use_fmincon!');
end

% mix fixed parameters in again
p.x(options.fixed == 0) = x;
x = p.x;

% restore warning settings
warning(s);

% functions ---------------------------------------------------------------

    function f = L(x)
        % Computes the neg. log. likelihood for given parameters x where
        % fixed parameters can be included. The model function, the noise
        % likelihood model and the likelihood penalty are specified in
        % variable p. Built-in functions can be used (see fit_func.m)

        % computes the likelihood of an n-D function with fixed and variable
        % parameters and a dataset (both given in the global variable p)

        % first mix fixed and variable parameters
        h = x;
        x = p.x;
        x(p.fix == 0) = h;

        % calculate the model function (external/built-in)
        y = fit_func(p.fun, x, p.xdata);

        % calculate the likelihood between model and data and add penalty
        f = Likelihood(p.lik, p.ydata, y);
        if isa(p.pen_fun, 'function_handle')
            f = f + p.pen_fun(x);
        end
    end

    function L = Likelihood(which, img, mod)
        % Computes neg. log. likelihood between image (img) and model (mod)
        % for the given noise model (which)
        switch which
            case 'Gaussian'
                % Gaussian negative logarithmic likelihood is proportional
                % to sum((img-mod).^2)
                L = (img - mod).^2;
            case 'Poisson'
                % Poisson negative logarithmic likelihood is defined as:
                % L = sum(mod - img * ln(mod) + ln(img!)) with the Poisson
                % probability P(v, n) = exp(-v)*v^n/n! of obtaining n(>=0)
                % events if the mean number of events is v
                % img must be a nonnegative integer array and mod must be
                % strictly positive
                % we set mod to greater than 1E-6 to avoid zeros.
                mod = max(1E-6, mod);
                L = mod - img .* log(mod);
                % however the additional summand ln(img!) is not dependent
                % on the parameters and therefore does not need to be
                % included in the minimization
            otherwise
                error('Likelihood type unsupported!');
        end
        L = sum(L(:));
        % finally sum up
    end

    function b = isnotfieldorisempty(s, fieldname)
        % Gives true if fieldname is not a field of struct s or
        % if s.fieldname is empty

        b = 1;
        if isfield(s, fieldname)
            if ~isempty(getfield(s, fieldname))
                b = 0;
            end
        end
    end

end

Error using ==> fit_nl_ex at 71
Not enough arguments given!

Error in ==> fit_nl_ex at 71
    error('Not enough arguments given!');

Published with MATLAB® 7.6