function [] = bootstrap2groupsMEANS (ctrl, exp)
%
% function [] = bootstrap2groupsMEANS (ctrl, exp)
%
% 'ctrl' and 'exp' are assumed to be column vectors, this is
% necessary for creation of 'box'
%
% This function will use bootstrap to compute the probability
% that data in 'ctrl' and 'exp' are from the same distribution,
% using the group means as descriptors, and the absolute value
% of the difference between means as the test statistic 
% (similar to Student's t-test)

% --- prepare for resampling ---

ctrl_d = mean(ctrl);                % compute descriptors                 
exp_d = mean(exp);                  % 'ctrl_d' and 'exp_d'

test_stat = abs(ctrl_d - exp_d)     % compute test statistic 'test_stat'   

box = [ctrl; exp];                  % combine datasets into 'box'
                                    % could use concat()
 
pseudo_stat_dist = zeros(1,10000);  % create vector 'pseudo_stat_dist'  
                                    % to store results of resampling

% --- start resampling ---

for trials = 1:10000                        % resample 10000 times
    
    pseudo_ctrl = ...                       % create pseudo groups to be
        sample(length(ctrl), box);          % same size as actual groups
    pseudo_exp = ...
        sample(length(exp), box);
    
    pseudo_ctrl_d = ...                     % compute pseudo group 
        mean(pseudo_ctrl);                  % descriptors
    pseudo_exp_d = ...
        mean(pseudo_exp);
    
    pseudo_stat = ...                       % compute pseudo test stat
        abs(pseudo_ctrl_d - pseudo_exp_d);     
    pseudo_stat_dist(trials) = ...          % store pseudo test stat to 
        pseudo_stat;                        % build null distribution
    
end

% --- end resampling ---  

bigger = count(pseudo_stat_dist > test_stat); % count pseudo test stats 
                                              % bigger than actual stat

pval = bigger / length(pseudo_stat_dist)      % divide by total number
                                              % of pseudo test stats to 
                                              % compute p-value
                                              % could use proportion()

hist(pseudo_stat_dist)                        % show histogram of pseudo
    title('Pseudo test stat distribution')    % test stats
    xlabel('Pseudo test stat')
    ylabel('Frequency')