bootstrapFeedbackStrength.m 2.8 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 78 79 80 81 82 83 84 85 86 87 88
function strength_values = bootstrapFeedbackStrength(data1,data2,sample_size,n_bootstraps,analysis_type)

% Bootstrap two time series and produce a historgram of
% slope parameters of a linear regression.
%
% USAGE: strength_values = bootstrapRegression(data1,data2,sample_size,n_bootstraps)
%
%     OUTPUT
%         strength_values: bootstrapped slope values
%
%     INPUT
%         data1:         Time series 1, "x-values" or predictors for the regression model
%         data2:         Time series 2, "y-values" or predictands of the regression model
%         sample_size:   How many samples should be drawn from the
%                        full data?
%         n_bootstraps:  How often should the bootstrapping be repeated?
%         analysis_type: Which analysis method should be used to diagnose the
%                        feedback strength? acc (ACC), reg (regular linear regression,
%                        that regresses >data2< on >data1<), or rreg (robust
%                        regression that regresses >data2< on >data1<)

% ----------- %
% PREPARATION %
% ----------- %

    % 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

% --------- %
% BOOTSTRAP %
% --------- %

    % total pool of data
    n_data       = length(data1);
    all_indices  = 1:n_data;

    % storage
    strength_values = zeros(1,n_bootstraps);

    % actual bootstrapping
    for c_strap = 1:n_bootstraps

        % randomly draw samples
        random_sample = datasample(all_indices,sample_size,'Replace',false);

        sample_data1 = data1(random_sample);
        sample_data2 = data2(random_sample);

        % calculate feedback strengths
        switch analysis_id

            % correlation
            case 1

                % acc
                current_feedback = corrcoef(sample_data1,sample_data2);
                current_feedback = current_feedback(1,2);

            % normal linear regression
            case 2

                regression_parameters = polyfit(sample_data1,sample_data2,1);
                current_feedback = regression_parameters(1);

            % robust regression
            case 3

                rregression_parameters = robustfit(sample_data1,sample_data2);
                current_feedback = rregression_parameters(2);  % yes, the order is different from polyfit ...

        % switch: calculate feedback strength
        end

        % assign regression value
        strength_values(c_strap) = current_feedback;

    % for: repeat bootstrapping
    end