nanmean gives unexpected results

I am calculating nanmean using two different techniques. I have a matrix avg_terra_aqua1 of size 28*360000. In the first method, I am calculating nanmean of the matrix along its first dimension as follow:

sp_aod1 = nanmean(avg_terra_aqua1, 1);
sp_aod1_modis = reshape(sp_aod1, [400, 900]);

Using the second method, I am calculating nanmean in a slightly more indirect way as follows:

nan_mat = isnan(avg_terra_aqua1);
count = squeeze(sum(nan_mat, 1));
  for j = 1:360000;
         if(count(1, j) < 29);                     
             sp_aod(1, j)  = nanmean(avg_terra_aqua1(1, j), 1);
         else
             sp_aod(1, j) = NaN;
         end;
  end;
  sp_aod_modis = reshape(sp_aod, [400, 900]);

I expect the results of sp_aod1_modis and sp_aod_modis to be similar, because in the second method, all the count values are below 29 (it is a dataset for the month of february which has 28 days). However I am getting very different results -- I have attached the results of two methods. First figure is the result of second method and second figure is the result of the first method. Can someone explain why this is happening?

Answer 1

[Stephen23]

"I am calculating nanmean using two different techniques"

Nope, not at all. Your code does not calculate anything at all like a mean, and it makes very little sense. Take a look at this line:

sp_aod(1, j)  = nanmean(avg_terra_aqua1(1, j), 1);

You use indexing to extract just one element from the first row of avg_terra_aqua1, and then take the nanmean of that one value.... which gives exactly the same value. What do you expect the average of one value to give you?

"I expect the results of sp_aod1_modis and sp_aod_modis to be similar"

I don't expect that. Your code simply stores the first row of the input data in a new variable sp_aod, so I expect that the output will just be the first row. Lets have a look with the help of some fake data:

avg_terra_aqua1 = randi(9,7,5);
avg_terra_aqua1(:,1) = NaN;
avg_terra_aqua1(3,:) = NaN

sp_aod1 = nanmean(avg_terra_aqua1, 1);

nan_mat = isnan(avg_terra_aqua1);
count = squeeze(sum(nan_mat, 1));
for j = 1:size(avg_terra_aqua1,2);
    if(count(1, j) < 29);
       sp_aod(1, j)  = nanmean(avg_terra_aqua1(1, j), 1);
    else
       sp_aod(1, j) = NaN;
    end
end

sp_aod1
sp_aod

displays this:

avg_terra_aqua1 =
   NaN     9     1     4     5
   NaN     9     2     8     6
   NaN   NaN   NaN   NaN   NaN
   NaN     6     6     8     5
   NaN     6     9     4     9
   NaN     6     2     9     6
   NaN     3     5     8     9
sp_aod1 =
     NaN   6.5000   4.1667   6.8333   6.6667
sp_aod =
   NaN     9     1     4     5
>>

As I expected, your loop simply returns the first row of the input matrix, without anything like a mean calculation. nanmean works correctly.

Summary: Your mistake is very simple, and is commonly made by beginners: you did not actually look at what the code is really doing. Beginners often come here after writing many lines of code and claim that it does "such and such". Nope, it doesn't. Those beginners made some mistake (that is oaky, we all make mistakes), but are stuck with thinking that just because they believe/want/claim/intend/... that their code really really should do "such and such", that therefore it must be doing that. The simplest way of checking if it really does that is to actually check that it does whatever it needs to do: check every line of code, make sure that it does what you need it to do, confirm every calculation by hand, run a small test case. Only when a line is correct should you move on to writing the next line.

[Sagar]

Thank you Stephen. That was my mistake. The code should have been:

sp_aod(1, j) = nanmean(avg_terra_aqua1(:, j), 1);

Now both the results look the same. I will leave this question undeleted.