Вычислить скользящее среднее данных рекурсивно

Question

У меня есть две двумерные матрицы A и B, где строки указывают на испытания, а столбцы указывают образцы, собранные в ходе испытания.

Я нахожусь в сценарии, где A доступно, но B собирается в режиме реального времени. Я хочу рассчитать текущее среднее значение {A и доступные данные для B}, так как B отбирается. Я думал, что смогу выполнить sh, рассчитав средневзвешенное значение A и B и обновив весовые коэффициенты по мере сбора проб и образцов для B. В частности, я подумал, что могу обновить веса и рекурсивно использовать значения, которые я уже сохранил из предыдущей итерации. Ниже приведен мой код и график вывода:

close all;
clear all;

%define the sizes of the matrices -- exact numbers aren't important for illustration
n1 = 5;
n2 = 10;
n3 = 12;

%define a matrix that will act as the history of data already collected
A = randi(10,[n2,n1]);
A_avg = mean(A,1); %averaged across n2 trials to get n1 values

%current acts as "incoming" data
B = randi(10,[n3,n1]); %n3 trials, n1 samples per trial

%preallocate matrices for final solutions
correct_means = zeros(n3,n1);
estimated_means = zeros(n3,n1);

for k1=1:size(B,1) %loop through trials
    %get running average in the case where we already have all samples
    correct_means(k1,:) = mean([A;B([1:k1],:)],1);
    for k2=1:size(B,2) %k2 should loop through samples
        %calculate averages as samples are incoming recursively (weighted averaging)
        if k1>1
            estimated_means(k1,k2) = (n2 / (n2+k1)) * A_avg(k2)...
                + ((k1-1)/(n2+k1)) * estimated_means(k1-1,k2) + (1/(n2+k1)) * B(k1,k2);
        elseif k1==1
            estimated_means(k1,k2) = (n2 / (n2+k1)) * A_avg(k2)...
                + ((k1-1)/(n2+k1)) * estimated_means(k1,k2) + (1/(n2+k1)) * B(k1,k2);       
        end
%         if k1==2, keyboard; end
    end
end

%plot the results
figure; hold on;
plot(nan, 'b', 'displayname', 'correct solution');
plot(nan, 'k--', 'displayname', 'my solution');
leg_tmp = legend('show');
set(leg_tmp,'Location','Best');

plot(correct_means, 'b', 'displayname', 'correct solution');
plot(estimated_means, 'k--', 'displayname', 'my solution');

ylabel('running averages');
xlabel('samples');

На прилагаемом графике изложены мои попытки решения (черный) и то, во что я верю быть правильным ответом (синий). Обратите внимание, что я строю среднее значение только после получения всех выборок для всех испытаний, но я сохраняю скользящее среднее по мере сбора данных. Как вы видите, мои ответы кажутся немного странными.

Моя идея заключалась в том, что A должно быть обновлено в соответствии с долей испытаний, использованных для определения его среднего по отношению к общему количеству испытаний по мере сбора B. Аналогично, вес для текущей выборки B просто равен 1, деленному на общее количество текущих испытаний в итерации, и предыдущие выборки B рекурсивно вызываются и взвешиваются соответственно. Эти веса в сумме равны 1 и имеют для меня смысл, поэтому мне трудно понять, где я все испортил.

Кто-нибудь может увидеть, где я все испортил?

Answer 1

Следует учесть, что чем длиннее код, тем больше ошибок он накапливает.

Заранее извиняюсь, что мне было проще переписать бизнес-логику c, чем найти ошибку в вашем коде - так что, хотя результат может быть не совсем то, что вы искали, он дает исправление. Я надеюсь, что вы все еще можете найти это полезным.

Пожалуйста, посмотрите на немного упрощенную версию, которая дает правильные результаты:

function q60180320
rng(60180320); % For reproducibility

% define the sizes of the matrices
n1 = 5;
n2 = 10;
n3 = 12;

% define a matrix that will act as the history of data already collected
A = randi(10,[n2,n1]);
A_avg = mean(A,1); %averaged across n2 trials to get n1 values

% current acts as "incoming" data
B = randi(10,[n3,n1]); %n3 trials, n1 samples per trial
correct_means = cumsum([A;B],1)./(1:n2+n3).'; correct_means = correct_means(n2+1:end,:);

% preallocate matrices
estimated_means = zeros(n3+1,n1);
estimated_means(1,:) = A_avg; % trick to avoid an if-clause inside the loop

for k1 = 1:size(B,1) % Loop through trials
  %% Compute weights:
  totalRows = n2 + k1;
  W_old = (totalRows - 1)./totalRows;
  W_new = 1/totalRows;

  %% Incoming measurement (assuming an entire row of B becomes available at a time)
  newB = B(k1,:); 

  %% Compute running average
  estimated_means(k1+1,:) = W_old * estimated_means(k1,:) + W_new * newB;
end
estimated_means = estimated_means(2:end,:); % remove the first row;

% plot the results
figure; hold on;
plot(nan, 'k', 'displayname', 'correct solution');
plot(nan, 'w--', 'displayname', 'my solution');
leg_tmp = legend('show');
set(leg_tmp,'Location','EastOutside');

plot(correct_means, 'k', 'displayname', 'correct solution', 'LineWidth', 2);
plot(estimated_means, 'w--', 'displayname', 'my solution');

ylabel('running averages');
xlabel('samples');