diff --git a/inst/boot.m b/inst/boot.m
index e0301101..e178a6fb 100755
--- a/inst/boot.m
+++ b/inst/boot.m
@@ -142,6 +142,7 @@
     % Assign weights (counts) for uniform sampling
     c = ones (n, 1) * nboot; 
   end
+  N = sum (c);
 
   % Perform balanced sampling
   r = 0;
@@ -156,15 +157,14 @@
       end
     end
     for i = 1:n
-      d = c;  
+      d = c;
       if (u)
-        d(r) = 0;
-      end
-      if (~ sum (d))
-        d = c;
+        if (N > d(r))
+          d(r) = 0;
+        end
       end
       d = cumsum (d);
-      j = sum (R(i) >= d ./ d(n)) + 1;
+      j = sum (R(i) * d(n) >= d) + 1;
       if (isvec) 
         bootsam (i, b) = x(j);
       else
@@ -172,11 +172,11 @@
       end
       if (nboot > 1)
         c(j) = c(j) - 1;
+        N = N - 1;
       end
     end
   end
 
-
 %!demo
 %!
 %! % N as input; balanced bootstrap resampling with replacement
diff --git a/inst/bootclust.m b/inst/bootclust.m
index f7a7bda2..ad2736e9 100755
--- a/inst/bootclust.m
+++ b/inst/bootclust.m
@@ -521,21 +521,31 @@
         warning ('off', 'all');
       end
       try
-        jackfun = @(i) bootfun (vertcat (x{1 : nx ~= i, :}));
-        if (ncpus > 1)  
-          % PARALLEL evaluation of bootfun on each jackknife resample 
-          if (ISOCTAVE)
-            % OCTAVE
-            T = cell2mat (pararrayfun (ncpus, jackfun, 1 : nx, ...
-                                       'UniformOutput', false));
+        if (VECTORIZED)
+          % Leave-one-out DATA resampling followed by vectorized function
+          % evaluations
+          T = bootfun (reshape (cell2mat (arrayfun (...
+                       @(i) vertcat (x{1 : nx ~= i, :}), (1 : nx)', ...
+                       'UniformOutput', false)), n - n / nx, []));
+        else
+          % Leave-one-out DATA resampling followed by looped function
+          % evaluations (if bootfun is not vectorized)
+          jackfun = @(i) bootfun (vertcat (x{1 : nx ~= i, :}));
+          if (ncpus > 1)  
+            % PARALLEL evaluation of bootfun on each jackknife resample 
+            if (ISOCTAVE)
+              % OCTAVE
+              T = cell2mat (pararrayfun (ncpus, jackfun, 1 : nx, ...
+                                         'UniformOutput', false));
+            else
+              % MATLAB
+              T = zeros (m, nx);
+              parfor i = 1 : nx; T(:, i) = feval (jackfun, i); end
+            end
           else
-            % MATLAB
-            T = zeros (m, nx);
-            parfor i = 1 : nx; T(:, i) = feval (jackfun, i); end
+            % SERIAL evaluation of bootfun on each jackknife resample
+            T = cell2mat (arrayfun (jackfun, 1 : nx, 'UniformOutput', false));
           end
-        else
-          % SERIAL evaluation of bootfun on each jackknife resample
-          T = cell2mat (arrayfun (jackfun, 1 : nx, 'UniformOutput', false));
         end
         % Calculate empirical influence function
         U = (nx - 1) * bsxfun (@minus, mean (T, 2), T);
diff --git a/inst/bootknife.m b/inst/bootknife.m
index 898dea02..6acf96d7 100755
--- a/inst/bootknife.m
+++ b/inst/bootknife.m
@@ -718,21 +718,31 @@
           % Attempt to form bias-corrected and accelerated (BCa) bootstrap CIs
           % Use the Jackknife to calculate the acceleration constant (a)
           try
-            jackfun = @(i) bootfun (x(1 : n ~= i, :));
-            if (ncpus > 1)  
-              % PARALLEL evaluation of bootfun on each jackknife resample 
-              if (ISOCTAVE)
-                % OCTAVE
-                T = cell2mat (pararrayfun (ncpus, jackfun, 1 : n, ...
-                                           'UniformOutput', false));
+            if (vectorized)
+              % Leave-one-out DATA resampling followed by vectorized function
+              % evaluations
+              T = bootfun (reshape (cell2mat (arrayfun (...
+                           @(i) x(1 : n ~= i, :), (1 : n)', ...
+                           'UniformOutput', false)), n - 1, []));
+            else
+              % Leave-one-out DATA resampling followed by looped function
+              % evaluations (if bootfun is not vectorized)
+              jackfun = @(i) bootfun (x(1 : n ~= i, :));
+              if (ncpus > 1)  
+                % PARALLEL evaluation of bootfun on each jackknife resample 
+                if (ISOCTAVE)
+                  % OCTAVE
+                  T = cell2mat (pararrayfun (ncpus, jackfun, 1 : n, ...
+                                             'UniformOutput', false));
+                else
+                  % MATLAB
+                  T = zeros (m, n);
+                  parfor i = 1 : n; T(:, i) = feval (jackfun, i); end
+                end
               else
-                % MATLAB
-                T = zeros (m, n);
-                parfor i = 1 : n; T(:, i) = feval (jackfun, i); end
+                % SERIAL evaluation of bootfun on each jackknife resample
+                T = cell2mat (arrayfun (jackfun, 1 : n, 'UniformOutput', false));
               end
-            else
-              % SERIAL evaluation of bootfun on each jackknife resample
-              T = cell2mat (arrayfun (jackfun, 1 : n, 'UniformOutput', false));
             end
             % Calculate empirical influence function
             if (~ isempty (strata))