Why do nested means give different answers depending on which dimension is averaged first

Consider that I have a square matrix where I want the average of all elements. Also consider that because I'm awkward I want to do this using mean() rather than mean2(). To do this I could take the average of rows then the average of that. Or I could do the same for columns. I've tried both here using the simple code below and the offset is always very small (e-16) but always non-zero. Am I missing something fundamental?

array = rand(100);

mean_1 = mean(mean(array,1));

mean_2 = mean(mean(array,2));

offset = mean_1 - mean_2;

[jonas]

Its probably because the first average gives an vector of floating point number which are not exact

Answer 1

[Adam]

It is just rounding errors. When I tried your code the offsets were of the order of 1e-16, which is close enough to 0 as to make no difference and is within the expected level of error for doing calculations in a different order, with rounding of the intermediate results based on how numbers are represented in floating point data types (many numbers are impossible to store to 100% accuracy as a float).

You can find many sources on this type of thing if you do an internet search for floating point data errors or similar phrases.

As an aside

mean( array(:) )

allows you to take the mean of the whole matrix in one go without mean2

[Barry Crouch]

Thanks, for this sanity now restored.