Цикл по элементам n
в x
и в цикле, вычисляющем дисперсию всех элементов до i
с использованием var(x(1:i))
, составляет алгоритм O (n 2 ). Это дорого по своей природе.
Выборочная дисперсия (которую вычисляет var
) определяется как sum((x-mean(x)).^2) / (n-1)
, с n = length(x)
. Это может быть переписано как (sum(x.^2) - sum(x).^2 / n) / (n-1)
. Эта формула позволяет нам накапливать sum(x)
и sum(x.^2)
в одном цикле, а затем вычислять дисперсию. Это также позволяет нам вычислить кумулятивную дисперсию в O (n).
Для вектора x
, у нас будет следующий цикл:
x = randn(100,1); % some data
v = zeros(size(x)); % cumulative variance
s = x(1); % running sum of x
s2 = x(1).^2; % running sum of square of x
for ii = 2:numel(x) % loop starts at 2, for ii=1 we cannot compute variance
s = s + x(ii);
s2 = s2 + x(ii).^2;
v(ii) = (s2 - s.^2 / ii) / (ii-1);
end
Мы можем избежать явногоцикл с использованием cumsum
:
s = cumsum(x);
s2 = cumsum(x.^2);
n = (1:numel(x)).';
v = (s2 - s.^2 ./ n) ./ (n-1); % v(1) will be NaN, rather than 0 as in the first version
v(1) = 0; % so we set it to 0 explicitly here
Код в OP вычисляет кумулятивную дисперсию для каждого столбца матрицы. Приведенный выше код может быть тривиально адаптирован для того же:
s = cumsum(VItimeseries,1); % cumulative sum explicitly along columns
s2 = cumsum(VItimeseries.^2,1);
n = (1:size(VItimeseries,1)).'; % use number of rows, rather than `numel`.
v = (s2 - s.^2 ./ n) ./ (n-1);
v(1,:) = 0; % fill first row with zeros, not just first element