bootstrapFeedbackStrengthDiff.m 3.55 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 89 90 91 92 93 94 95 96 97 98 99
function strength_diffs = bootstrapFeedbackStrengthDiff(data1,data2,sample_size,n_bootstraps,analysis_type)

% Bootstrap tuples of two time series two time series and produce a historgram of
% of differences of the feedback strength diagnosed from either ACC, linear regression
% or robust regression. The differences are calculated
% as follows: Randly sample m tuples from the time series and calculate the slope
% parameter; then use the remaining n-m tuples, again calculate the slope parameter,
% and save the difference between the two. n is total length of the time dataset,
% and m is equivalent with the parameter sample_size.
%
% USAGE: strength_diffs = bootstrapRegression(data1,data2,sample_size,n_bootstraps)
%
%     OUTPUT
%         strength_diffs: bootstrapped slope difference 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_diffs = zeros(1,n_bootstraps);

    % actual bootstrapping
    for c_strap = 1:n_bootstraps

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

        % opposite sample
        opposite_sample = setdiff(all_indices,random_sample);

        % calculate feedback strength differences
        switch analysis_id

            % correlation
            case 1

                % acc
                random_feedback   = corrcoef(data1(random_sample),    data2(random_sample));
                opposite_feedback = corrcoef(data1(opposite_feedback),data2(opposite_feedback));

                feedback_difference = random_feedback(1,2)-opposite_feedback(1,2);

            % normal linear regression
            case 2

                random_pp   = polyfit(data1(random_sample),data2(random_sample),1);
                opposite_pp = polyfit(data1(opposite_sample),data2(opposite_sample),1);

                feedback_difference = random_pp(1)-opposite_pp(1);

            % robust regression
            case 3

                random_rreg   = robustfit(data1(random_sample),data2(random_sample));
                opposite_rreg = robustfit(data1(opposite_sample),data2(opposite_sample));

                feedback_difference = random_rreg(2)-opposite_rreg(2);

        % switch: calculate feedback strength
        end

        % assign regression value
        strength_diffs(c_strap) = feedback_difference;

    % for: repeat bootstrapping
    end