function ret = pairwise_mean (x, dim) siz = size (x); if (nargin < 2) dim = find (siz > 1,1); if isequal (dim,[]) dim = 1; endif endif n = siz(dim); nn = n; siz(dim) = 1; L = prod (size (x)(1:dim-1)); inds = ((1:L)'+L*n*(0:prod (siz)/L-1))(:); k = 0; while n > 1 # x(inds+L*(0:n/2-1)) += x(inds+L*(ceil (n/2):n-1)); # x(inds+L*(0:(n-1)/2)) /= 2; x(inds+L*(0:n/2-1)) += (x(inds+L*(ceil (n/2):n-1)) - x(inds+L*(0:n/2-1))) / 2; x(inds+L*floor(n/2)) /= 2; n = ceil (n/2); ++k; end ret = reshape (x(inds),siz) * (2^k / nn); %!test %! x = randn (10, 11, 12); %! assert (pairwise_mean (x), mean (x), 2*eps) %!test %! x = randn (10, 11, 12); %! assert (pairwise_mean (x, 1), mean (x, 1), 2*eps) %!test %! x = randn (10, 11, 12); %! assert (pairwise_mean (x, 2), mean (x, 2), 2*eps) %!test %! x = randn (10, 11, 12); %! assert (pairwise_mean (x, 3), mean (x, 3), 2*eps) %!test %! x = [1.7e308, 1.5e308]; %! assert (pairwise_mean (x), 1.6e308, 2*eps) %!test %! rand ("state",3); x = rand (10, 1) * 1e-323; %! assert (pairwise_mean (x) > 0, true)