Why does bsxfun produce different result than brute force mean?

Question

KAE 2017-10-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/362786-why-does-bsxfun-produce-different-result-than-brute-force-mean

评论： KAE 2017-10-23

I would like to apply a 2D mask to a 3D matrix, then average the non-masked values. In the example below, the mean is calculated in two ways, (a) with bsxfun and (b) in a loop which sums values and divides by the count. The results aren't the same, though. Can you spot my error, which I suspect is in approach (b)?

nRow = 200; nCol = 100; nChannel = 10;
bigMatrix = rand(nRow, nCol, nChannel); % 3D matrix of values that we want to mask
mask = ones(nRow, nCol); % Make a 2D mask
mask(1:25, 26:50) = NaN; % Mask set to NaN where values should be excluded
%%(a) Use bsxfun to substitute NaNs for masked values in bigMatrix
% Since mask values are either NaN or 1, can mask by multiplying 
maskedByMult = bsxfun(@times, bigMatrix, cast(mask, 'like', bigMatrix));
% Average the non-nan values
meanByBsxfun = squeeze(mean(mean(maskedByMult, 'omitnan'), 'omitnan')); 
%%(b) Mask by loop
nGoodValues = 0; % Count of unmasked values used in later calculation of their mean
sumValues = zeros(nChannel, 1); % Sums will be added to this
for iRow = 1:length(nRow)
    for iCol = 1:length(nCol)
        if isfinite(mask(iRow, iCol)) % Check if this x,y location is masked out
            nGoodValues = nGoodValues + 1;
            sumValues = sumValues + squeeze(bigMatrix(iRow, iCol, :));
        end
    end
end
meanByLoop = sumValues/nGoodValues; % Get the mean
%%Show they differ
figure('Name', 'Compare two ways of getting the mean');
plot(1:10, meanByBsxfun); hold on;
plot(1:10, meanByLoop, '--');
legend('Using bsxfun', 'Using a loop');

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2017-10-23

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/362786-why-does-bsxfun-produce-different-result-than-brute-force-mean#answer_287290

编辑：David Goodmanson 2017-10-23

在 MATLAB Online 中打开

Hello KAE,

I believe that each method has a small problem. In the second method, since length(nRrow) = length(nCol) = 1, the for loops only execute one time. So you need

for iRow = 1:nRow % etc.

In the first method, denote maskedbyMult by A, and let N = the number of elements that are not nan, which in this case is 200*100 - 625. Then

mean(mean(A),'omitnan')   is not equal to    sum(sum(A),'omitnan')/N    which is what you want.

for example,

a = magic(3);  
a(1,2) = nan
a =
     8   NaN     6
     3     5     7
     4     9     2
mean(mean(a,'omitnan'))
ans = 5.6667
sum(sum(a,'omitnan'))/8
ans = 5.5000

After those changes both methods agree, and both agree with

bigMatMask = repmat(mask,1,1,10).*bigMatrix;
mean3 = squeeze(sum(sum(bigMatMask,2,'omitnan'),1)/(200*100-625))
% if you have a newer version of Matlab with explicit expansion,
% the first line can be    bigMatMask = mask.*bigMatrix;

3 个评论
显示 1更早的评论隐藏 1更早的评论

KAE 2017-10-23

编辑：KAE 2017-10-23

在 MATLAB Online 中打开

Thanks so much, since I have used mean(mean(x)) erroneously in lots of other code. By the way here are the time differences on my machine,

bxsfun (a) 0.003004 s
loop (b) 0.108191 s
repmat without implicit expansion since I am using R2015b (c) 0.003103 s

And here is the corrected full code if that helps anyone.

nRow = 200; nCol = 100; nChannel = 10;
bigMatrix = rand(nRow, nCol, nChannel); % 3D matrix of values that we want to mask
mask = ones(nRow, nCol); % Make a 2D mask
mask(1:25, 26:50) = NaN; % Mask set to NaN where values should be excluded
nIn = nRow*nCol - sum(isnan(mask(:))); % How many non-masked values are in matrix
%%(a) Use bsxfun to substitute NaNs for masked values in bigMatrix
tic;
% Since mask values are either NaN or 1, can mask by multiplying
maskedByMult = bsxfun(@times, bigMatrix, cast(mask, 'like', bigMatrix));
% Average the non-nan values
meanByBsxfun = squeeze(sum(sum(maskedByMult, 'omitnan'))/nIn);
toc;
%%(b) Mask by loop
tic;
nGoodValues = 0; % Count of unmasked values used in later calculation of their mean
sumValues = zeros(nChannel, 1); % Sums will be added to this
for iRow = 1:nRow
    for iCol = 1:nCol
        if isfinite(mask(iRow, iCol)) % Check if this x,y location is masked out
            nGoodValues = nGoodValues + 1;
            sumValues = sumValues + squeeze(bigMatrix(iRow, iCol, :));
        end
    end
end
meanByLoop = sumValues/nGoodValues; % Get the mean
toc;
%%(c) Repmat
tic;
bigMask = repmat(mask,1,1,nChannel).*bigMatrix;
meanByRepmat = squeeze(sum(sum(bigMask,2,'omitnan'),1)/nIn);
toc;
%%Show they are the same
figure('Name', 'Compare two ways of getting the mean');
plot(1:10, meanByBsxfun); hold on;
plot(1:10, meanByLoop, '--');
plot(1:10, meanByRepmat, ':');
legend('Using bsxfun', 'Using a loop', 'Using repmat');

David Goodmanson 2017-10-23

nice to see the code comparison.

mean(mean( )) of course does work out in lots of cases such as rectangular matrices without any nans.

By the way I meant to say implicit expansion, not explicit expansion.

KAE 2017-10-23

在 MATLAB Online 中打开

If someone wants to do another operation besides mean on a masked matrix, here is an example using median,

        % Reshape matrix into [spatial dimension x channel]
        temp = reshape(maskedByMult, [nRow*nCol nChannel]); 
        medianByBsxfun = median(temp, 'omitnan');

请先登录，再进行评论。