calculateFeedback.m 2.51 KB
Newer Older
Claas Faber's avatar
Claas Faber committed
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 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
function feedback_strength = calculateFeedback(forcing,response,analysis_type)

% Calculate the "feedback strength" (co-variability) between two variables
%
% SYNTAX:
%
%     feedback_strength = calculateFeedback(forcing,response,analysis_type)
%
%     OUTPUT:
%         feedback_strength: Scalar, "feedback strength" between the >forcing< and
%                            and >response< datasets, based on the chosen >analysis_type<.
%
%
%     INPUT:
%         forcing:         Vector, forcing dataset ("forcing" of the "feedback").
%         response:               Vector of the same size as >forcing<, "response" of the "feedback"
%         analysis_type:     Char, specifies how the "feedback" strength is diagnosed.
%                            Either of
%                            "acc":  Anomaly correlation, not suitable for very short datasets
%                            "reg":  reguar linear regression, where >response< is regressed on >forcing<
%                            "rreg": robust lineat regression, where >response< is regressed on >forcing<

% --------------- %
% CHECK ARGUMENTS %
% --------------- %

    % do the datasets have the same size?
    size_forcing  = size(forcing);
    size_response = size(response);

    same_size = sum(~( size_forcing == size_response ));

    if same_size ~= 0
        error('>forcing< and >response< datasets do not have the same size.');
    end

    % assign an ID to the analysis type
    % for use with switch
    if strcmp(analysis_type,'acc')
        analysis_id  = 1;
    elseif strcmp(analysis_type,'reg')
        analysis_id  = 2;
    elseif strcmp(analysis_type,'rreg')
        analysis_id  = 3;
    else
        error_string = sprintf('Analysis type >%s< not known.\nChoose either of acc, reg, or rreg.',analysis_type);
        error(error_string);
    end


% --------------------------- %
% CALCULATE FEEDBACK STRENGTH %
% --------------------------- %

    switch analysis_id

        % correlation
        case 1

            % acc
            feedback_strength = corrcoef(forcing,response);
            feedback_strength = feedback_strength(1,2);

        % normal linear regression
        case 2

            regression_parameters = polyfit(forcing,response,1);
            feedback_strength = regression_parameters(1);

        % robust regression
        case 3

            rregression_parameters = robustfit(forcing,response);
            feedback_strength = rregression_parameters(2);  % yes, the order is different from polyfit ...

    % switch: calculate feedback strength
    end