rri_boot_order.m

国外的一个PLStoolbox,主要用于处理图象,也可以用来回归,欢迎使用

function [boot_order, new_num_boot] ...
   = rri_boot_order(num_subj_lst, num_cond, num_boot, incl_seq, ...
	min_subj_per_group, is_boot_samples, boot_samples, ...
	new_num_boot, nopopup)
%
%  USAGE:  [boot_order, new_num_boot] = rri_boot_order(num_subj_lst, num_cond, num_boot [, incl_seq = 0])
%   
%  Generate bootstrap sample order for all conditions of the groups.  
%  Assume the data rows are stored in the order of subject in condition
%  in group.
%   
%  NOTE: the bootstrap order are applied for the subjects within each
%        group only, i.e. sampling with replacement for the subjects
%        within group.  In order to make sure bootstrap results are
%        reliable, there must have at least 3 subjects in any of the
%        groups.
%
%  Method:  
%     The resampling orders are performed for each group separately.
%     The orders are generated by using MATLAB's 'rand' function. 
%
%  Input:
%         num_subj_lst: a N elements vector, each element specifies the
%			number of subjects of one of groups (subj_group).
%         num_cond:	number of conditions in the experiment
%			(assume the st_datamat has 1 row for 
%			each condition of each group).
%         num_boot:	number of resampling to be performed.
%         incl_seq:	switch to control whether the boot_order includes
%			a sequential order (original un-resampled order).
%			0 (default) means that sequential order should not
%			be included.
%			1 means that sequential order should be included.
%
%  Output: 
%         boot_order:	an MxN matrix that stores the new order for the N 
%			permutations, where M is the total number of rows 
%			of the st_datamat.
%			(i.e M = sum(num_subj_lst) * num_cond)
%         new_num_boot:	revised number of bootstrap resampling
%

   if ~exist('nopopup','var')
      nopopup = 0;
   end

   if ~exist('incl_seq','var')
      incl_seq = 0;
   end

   total_subj = sum(num_subj_lst);
   num_group = length(num_subj_lst);
   total_rows = num_cond * total_subj;

   if exist('min_subj_per_group','var') & isempty(min_subj_per_group)
      new_num_boot = num_boot;
      boot_order = [];
      return;
   end


if ~exist('min_subj_per_group','var') | isempty(min_subj_per_group)

   new_num_boot = num_boot;
   percentage = [];
   done = 0;

   %  num_subj in one of group is less than 3
   %
   if (min(num_subj_lst) < 3)

%      if exist('progress_hdl','var') & ishandle(progress_hdl)
%         close(progress_hdl);
%      end

      if nopopup
         disp('Bootstrap analysis requires that each group must have at least 3 subjects');
      else
         msg = 'Bootstrap analysis requires that each group must have at least 3 subjects';
         uiwait(msgbox(msg,'ERROR','modal'));
      end
%      disp('At least 3 subjects must be chosen for bootstrap analysis.');

      boot_order = [];
      return;
   end;

   while ~done
      if nopopup
         disp(' ');disp('Please input a percentage number (between 30 and 70) that will represent a');
         disp('minimum number of different subjects in bootstrap versus total number of');
         msg = 'subjects you have. If you are not sure, please use 50 percent: ';
         tmp=input(msg, 's');disp(' ');

         if isempty(tmp)
%            if exist('progress_hdl','var') & ishandle(progress_hdl)
%               close(progress_hdl);
%            end

            boot_order = [];
            return;
         else
            tmp = {tmp};
         end
      else
         msg = 'Please input a percentage number (between 30 and 70) that will represent a minimum number of different subjects in bootstrap versus total number of subjects you have. If you are not sure, please accept the default value which is equal to 50 percent:';
         tmp = inputdlg({msg}, 'Percentage of minimum different subject number', 1, {'50'});

         if isempty(tmp{1})
%            if exist('progress_hdl','var') & ishandle(progress_hdl)
%               close(progress_hdl);
%            end

            boot_order = [];
            return;
         end
      end

      percentage = str2num(tmp{1});

      if isempty(percentage) | percentage < 30 | percentage > 70
         done = 0;
      else
         done = 1;
      end
   end

   min_subj_per_group = ceil(min(num_subj_lst)*percentage/100);
   max_subj_per_group = 8;
%   tmp_boot_order = [];
%   progress_hdl = rri_progress_ui('initialize');

   %  if num of subj is 8 or less, we can get from boot matrix
   %
   is_boot_samples = zeros(1,num_group);
   boot_samples = cell(1,num_group);
   rand('state',sum(100*clock));

   if (sum(num_subj_lst <= max_subj_per_group) == num_group)

      for g = 1:num_group

         num_subj = num_subj_lst(g);
%         diff_subj = min_subj_per_group;

         boot_sample2 = [];

         %  remove sequential order
         %
         if incl_seq
            max_diff_subj = num_subj;
         else
            max_diff_subj = num_subj - 1;
         end

         for diff_subj = min_subj_per_group:max_diff_subj
            boot_sample1 = rri_boot_samples(num_subj, diff_subj);
            boot_sample2 = [boot_sample2; boot_sample1];
         end

         num_boot_samples = size(boot_sample2, 1);

         while num_boot_samples < new_num_boot

            done = 0;
            while ~done
               if nopopup
                  disp(' ');disp([num2str(num_subj), ' subjects can only have ',num2str(num_boot_samples), ...
			' different bootstrap samples. Please reduce the']);
                  disp('number of bootstrap (you also have option to choose another percentage');
                  msg1 = 'number later): ';
                  tmp1=input(msg1, 's');disp(' ');

                  disp(' ');disp('If you fail to reduce the number of bootstrap, please reduce percentage');
                  disp('number now (between 30 and 70) that will represent a minimum number of');
                  disp('different subjects in bootstrap versus total number of subjects you have.');
                  msg2 = 'If you are not sure, please use 50 percent: ';
                  tmp2=input(msg2, 's');disp(' ');

                  if isempty(tmp1) | isempty(tmp2)
%                     if exist('progress_hdl','var') & ishandle(progress_hdl)
%                        close(progress_hdl);
%                     end

                     boot_order = [];
                     return;
                  else
                     tmp = {tmp1, tmp2};
                  end
               else
                  msg1 = [num2str(num_subj), ' subjects can only have ',num2str(num_boot_samples), ...
			' different bootstrap samples. Please reduce the number of bootstrap:'];
                  msg2 = 'Or, please reduce percentage number (between 30 and 70) that will represent a minimum number of different subjects in bootstrap versus total number of subjects you have. If you are not sure, please accept the default value which is equal to 50 percent:';
                  tmp = inputdlg({msg1, msg2}, 'Edit', 1, {num2str(new_num_boot), num2str(percentage)});

                  if isempty(tmp{1}) | isempty(tmp{2})
%                     if exist('progress_hdl','var') & ishandle(progress_hdl)
%                        close(progress_hdl);
%                     end

                     boot_order = [];
                     return;
                  end
               end

               new_num_boot = str2num(tmp{1});
               percentage = str2num(tmp{2});
               min_subj_per_group = ceil(min(num_subj_lst)*percentage/100);

               if isempty(percentage) | percentage < 30 | percentage > 70 | new_num_boot < 1
                  done = 0;
               else
                  done = 1;
               end
            end

            if isempty(tmp)
%               if exist('progress_hdl','var') & ishandle(progress_hdl)
%                  close(progress_hdl);
%               end

               if nopopup
                  disp('You need more subjects in order to run bootstrap again.');
               else
                  msg = 'You need more subjects in order to run bootstrap again.';
                  uiwait(msgbox(msg,'ERROR','modal'));
               end

               boot_order = [];
               return;
            end

            boot_sample2 = [];

            %  remove sequential order
            %
            if incl_seq
               max_diff_subj = num_subj;
            else
               max_diff_subj = num_subj - 1;
            end

            for diff_subj = min_subj_per_group:max_diff_subj
               boot_sample1 = rri_boot_samples(num_subj, diff_subj);
               boot_sample2 = [boot_sample2; boot_sample1];
            end

            num_boot_samples = size(boot_sample2, 1);

         end	% while

         if ~isempty(boot_sample2)
            boot_sample2 = boot_sample2(randperm(num_boot_samples),:);
            boot_samples{g} = boot_sample2;
            is_boot_samples(g) = 1;
         end

      end	% for
   end		% if

end


   %  determine tmp_boot_order, which is re-ordered subject index matrix
   %
   tmp_boot_order = zeros(total_subj, new_num_boot);
   cnt = 0;

   for p=1:new_num_boot,

      subj_order = cell(1,num_group);
      not_done = 1;

      while (not_done)

         start_subj = 1;

         for g = 1:num_group
            num_subj = num_subj_lst(g);

            %  reorder tasks for the current group.
            all_samples_are_same = 1;

            while (all_samples_are_same)

               if is_boot_samples(g)		% get from boot_samples

                  new_subj_order = boot_samples{g}(p,:);
                  all_samples_are_same = 0;

                  not_done = 0;

               else

                  not_done = 1;

                  new_subj_order = floor(rand(1,num_subj)*num_subj) + 1;

                  % first_subj = new_subj_order(1);
                  % if ~isequal(new_subj_order,repmat(first_subj,1,num_subj))

                  test=length(unique(new_subj_order));

                  %% check to make sure there are more than n/2 people
                  %%
                  %if (test >= num_subj/2)

                  % check to make sure there are at lease min_subj_per_group people
                  %
                  if (test >= min_subj_per_group)
                     all_samples_are_same = 0;
                  end;
               end
            end;

            subj_order{g} = new_subj_order + start_subj - 1; 
            start_subj = start_subj + num_subj;

         end;

         if ~all(is_boot_samples)

            %  make sure the order is not a repeated one
            not_done = 0;
            for i=1:p-1,
               if isequal(squeeze(tmp_boot_order(:,i)),subj_order)
                  not_done = 1;
                  break;
               end;
            end;

           %  treat sequential order as duplicated one
           %
           if ~incl_seq & isequal([1:total_rows], [subj_order{:}])
              not_done = 1;
           end

            %  discard if all elements are the same
            cnt = cnt+1;
            if (cnt > 500),
               not_done = 0;
               disp('ERROR:  Duplicated bootstrap orders are used!');
            end;

         end;		% if ~all(is_boot_samples)

      end;		% while (not_done)

      tmp_boot_order(:,p) = [subj_order{:}]';

   end;  % for new_num_boot

  %  construct the resampling order matrix for bootstrap
  %
%  row_idx = reshape([1:total_rows],num_cond,total_subj);

  first = 1;
  last = 0;
  row_idx = [];
  for g=1:length(num_subj_lst)
     last = last + num_cond*num_subj_lst(g);
     tmp = reshape([first:last],num_subj_lst(g),num_cond);
     row_idx = [row_idx, tmp'];
     first = last+1;
  end

  boot_order = zeros(num_cond,total_subj, new_num_boot);

  for p=1:new_num_boot
     boot_order(:,:,p) = row_idx(:, tmp_boot_order(:,p));
  end

  b_order = [];
  for g=1:length(num_subj_lst)

     one_group = [];
     for p=1:new_num_boot
        tmp = ...
           boot_order(:,[sum(num_subj_lst(1:(g-1)))+1:sum(num_subj_lst(1:g))],p);
        tmp = reshape(tmp', [num_cond*num_subj_lst(g),1]);
        one_group = [one_group, tmp];
     end

     b_order = [b_order; one_group];

  end

  boot_order = b_order;

  return;