function ret = iter_mean (x, dim)
  if (nargin < 2)
    if (isvector (x))
      dim = columns (x);
    else
      dim = 1;
    endif
  endif
  l = length (size (x));
  S.subs = cell(1,l);
  S.subs(1:l) = ":";
  S.subs{dim} = 1;
  S.type = "()";
  ret = subsref (x, S);
  n = size (x, dim);
  for i = 2:n
    S.subs{dim} = i;
    ret += (subsref (x, S) -ret) / i;
  end
%!test
%! x = randn (10, 11, 12);
%! assert (iter_mean (x), mean (x), 2*eps)
%!test
%! x = randn (10, 11, 12);
%! assert (iter_mean (x, 1), mean (x, 1), 2*eps)
%!test
%! x = randn (10, 11, 12);
%! assert (iter_mean (x, 2), mean (x, 2), 2*eps)
%!test
%! x = randn (10, 11, 12);
%! assert (iter_mean (x, 3), mean (x, 3), 2*eps)
%!test
%! x = [1.7e308, 1.5e308];
%! assert (iter_mean (x), 1.6e308, 2*eps)
%!test
%! x = rand (10, 1) * 1e-323;
%! assert (iter_mean (x) > 0, true)
