Flip-Angle-Design-Toolbox/parameter_estimation.m at master · LarsonLab/Flip-Angle-Design-Toolbox · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
function [parameters_of_interest_est, nuisance_parameters_est] ...
    = parameter_estimation(y, model, goodness_of_fit_criterion, thetas)
    % FUNCTION PARAMETER_ESTIMATION computes an estimate of model parameters from data y
    %   Choice of goodness_of_fit_criterion:
    %      * 'least-squares'
    %      * 'maximum-likelihood'
    %
    % Flip Angle Design Toolbox
    % John Maidens (maidens@eecs.berkeley.edu)
    % June 2014

    display('===== Estimating parameter values =====')

    if strcmp(goodness_of_fit_criterion, 'least-squares')

        % perform least-squares fit
        options = optimset('MaxFunEvals', 5000, 'MaxIter', 5000);
        % obj = @(p) least_squares_objective(p, y, thetas*model.flip_angle_input_matrix', model);
        obj = @(p) least_squares_objective(p, y, thetas, model);
        p_opt = lsqnonlin(obj, ...
            [model.parameters_of_interest_nominal_values, ...
            model.nuisance_parameters_nominal_values], [], [], options);

        % unpack estimates
        l = length(model.parameters_of_interest_nominal_values);
        parameters_of_interest_est = p_opt(1:l);
        nuisance_parameters_est = p_opt(l+1:end);


    elseif strcmp(goodness_of_fit_criterion, 'maximum-likelihood')
        if strcmp(model.noise_type, 'Rician')

            % perform maximum-likelihood fit
            options = optimoptions('fminunc', 'MaxFunEvals', 5000, ...
                'MaxIter', 5000, 'Algorithm', 'quasi-newton');
            obj = @(p) negative_log_likelihood_rician(p, y, thetas, model);

            % initial point
            p_init =  [model.parameters_of_interest_nominal_values, ...
                model.nuisance_parameters_nominal_values];
            p_opt = fminunc(obj, p_init, options);

            % unpack estimates
            l = length(model.parameters_of_interest_nominal_values);
            parameters_of_interest_est = p_opt(1:l);
            nuisance_parameters_est = p_opt(l+1:end);

        else
            error(['Cannot do maximum likelihood fit for noise of type "', ...
                model.noise_type ,'"'])
        end


    else
        error(['Goodness of fit criterion "' , goodness_of_fit_criterion , ...
            '" is invalid'])
    end


end