CUMMEAN along elements of 3D matrix

I am trying to use cummean in order to take the cumulative mean by element in a (192x144x8766) matrix A.

I have tried cummean(A,3) and cummean(:,:,3) and get difference answers. Which one would be correct for taking the cumulative mean along the elements of a 3D matrix?

Answer 1

[Walter Roberson]

cummean(:,:,3) would be for accessing a variable along the third dimension, not for calling a function.

You should probably use cummean(A,3) . Except, that is, for the fact that cummean() is not a Mathworks function. There is a cummean() in the File Exchange. Or you can implement it yourself:

Acm3 = bsxfun(@divide, cumsum(A, 3), reshape(1:size(A,3), 1, 1, []) );

Or, if you have R2016b or later,

Acm3 = cumsum(A,3) ./ reshape(1:size(A,3), 1, 1, []);

[Vesp]

Hello Walter,

Thank you for your response. I am currently using cummean from the File Exchange for calculating the x and y components of IVT.

What I do is sort total ivt in descending order and then sort the x and y components like that of the sorted total ivt. When I do the cummean of ivtx and ivty then take their magntiude it does not equal the cummean of the sorted total ivt, which should be the case. Instead I get difference at each element for the final matrices of about 100. Would the error be this large if I am using cummean for matrices of AxBx8766 ?

[Walter Roberson]

What is IVT in this context? How is it calculated from x and y?

Is it correct that you have something like,

ivt = IVT(X, Y);

[sorted_ivt, sortidx] = sort(ivt, 3, 'descend');
ivt_cummean = cummean(sorted_ivt, 3)

ivtx = X(:, :, sortidx);
ivtx_cummean = cummean(ivtx, 3);

ivty = Y(:, :, sortidx);
ivty_cummean = cummean(ivty, 3);

ivt_from_cummeans = IVT(ivtx_cummean, ivty_cummean);

and you are looking at

difference = sorted_ivt - ivt_from_cummeans;

and finding that abs(difference) is approximately 100 for each location?

[Vesp]

In my case, IVT is integrated vapor transport and it is calculated the following way:

IVT=sqrt(x.^2+y.^2);

The process you have mentioned above is simular to mine. Please see a sample code for ivtx below (it includes extracting the top values defined by another matrix called ARcount):

totivt= sqrt((ivtx1).^2+ (ivty1).^2);
[sortedIVT,idx] =sort(totivt,3,'descend');

ivtx1(isnan(ivtx1))=0; %want to make sure NaN is not the maximum when sorting so make it zero
ivty1(isnan(ivty1))=0;
%need ivtx to be sorted in the order of total IVT and then take the

%sort ivtx like that of total ivt
[m,n,o]=size(ivtx1);
[Col,Row,~]=meshgrid(1:n,1:m,1:o);
lix = sub2ind([m,n,o],Row,Col,idx);
ivtxsort= ivtx1(lix);   %sort ivtx like total ivt was sorted 
     [x,y,z]=size(ivtxsort);
     t = ARcount > 0;
     C = zeros(x,y);
    [ii,jj] = ndgrid(1:x,1:y);
    C0= cummean(ivtxsort,3) %take the cumulative mean *%%%SOMETHING HERE GOES WRONG doesnt equal TOTIVT*
    C(t) = C0(sub2ind([x,y,z],ii(t),jj(t),ARcount(t))); %only apply the top N defined in historical 
    C(C==0)=NaN;
    Topivtx=C;

The same is done for ivt y component. For the sorted total ivt I take the cumulative mean like so:

      cumivt = cummean(sortedIVT,3);

To check I do the following:

    test=sqrt(C0.^2+C2.^2);
     isequal(cumivt,test) %result is 0 
     diff=cumivt-test; %the values here vary at each element from 0.0043 to 100+

Thank you for your time on this. I have been trying to figure this out for about a month.

[Image Analyst]

Just to check...what are you using cummean for? By chance, you don't want the mean of elements in a sliding/moving window do you? Because that would be movmean(). Or conv() or conv2() or convn().

[Vesp]

I am not sure what you mean by mean of elements in a sliding/moving window. I hadn't considered that as an option for what I would like this code to do.

I am using cummean as an alternative to selecting the Top nth values (defined by a AxB matrix), summing them and dividing by the number of N values defined by AxB matrix at each element in the 3D matrix. I figured if I did cummean then I could just use the AxB matrix to index the nth value specified for each element using this line:

C(t) = C0(sub2ind([x,y,z],ii(t),jj(t),ARcount(t)));

To make sure I am clear here is an example:

ARcount at (1x20)= 5
sortedivtx=(50x50x8766) 
cummean(sortedivtx)= 50x50x8766

but at 1x20 I would like the 5th cummean value defined by the ARcount matrix

[Walter Roberson]

The outline code I posted above turns out to have bugs in the construction of the sorted x and y data. The following is a working example

%create some data to test with
r = 3; c = 5; p = 7;
X = randi(50, r, c, p);
Y = randi(50, r, c, p);

IVT = @(x,y) x.^3 - x .* y.^2 - y + 3;  %create function to test with

proj_dim = 3;           %dimension to project against

cummean = @(A,dim) cumsum(A,dim) ./ permute(1:size(A,dim),circshift(1:ndims(A),dim-2));

ivt = IVT(X, Y);

[~, sortidx] = sort(ivt, proj_dim, 'descend');
[tr, tc, tp] = ndgrid(1:r, 1:c, 1:p);
idx_cell = {tr, tc, tp};
idx_cell{proj_dim} = sortidx;
sortind = sub2ind([r,c,p], idx_cell{:});

sorted_ivt = ivt(sortind);
ivt_cummean = cummean(sorted_ivt, proj_dim);

ivtx = X(sortind);
ivtx_cummean = cummean(ivtx, proj_dim);

ivty = Y(sortind);
ivty_cummean = cummean(ivty, proj_dim);

ivt_from_cummeans = IVT(ivtx_cummean, ivty_cummean);

difference = sorted_ivt - ivt_from_cummeans;

[~, maxind] = max( abs(difference(:)) );

fprintf('Maximum difference is %g\n', difference(maxind))
fprintf('Which occurs at X = %g, Y = %g original coordinates\n', ivtx(maxind), ivty(maxind) );
fprintf('Which is at X = %g, Y = %g cummean coordinates\n', ivtx_cummean(maxind), ivty_cummean(maxind) );

You are wondering what the cause of the difference is. It is because the change in function value for this sample IVT function is not linear with respect to changes in the coordinates.

Consider for example the function 1./x.^2 . Suppose you sample it at x = -1 and x = +1; you get an output of 1 for each of the two. cummean of [1 1] is [1 1] -- so you are projecting that if you sample at cummean([-1 1]) that you will get [1 1] as output. But cummean([-1 1]) is [-1 0], the outputs for which are [1 1/0] so the projection is infinitely wrong because of the nonlinear nature of the function.

Meanwhile, please go back and re-check how you sorted your 3D variables into consistent order. I show a method above that builds the appropriate linear indexing array.

The above code would have to be touched up a bit if you went to 4 or more dimensions, but not much really.

[Vesp]

Wow, it's incredible how much of a maximum difference there can be. I ran the code you provided above with my data and I got Maximum of 1305.47 and then checked the maximum difference with the original code I was using and I got a maximum difference of 1323.4.

Is there any way to reduce this difference? I have tried this by cumsum at each elements in the 3D matrix, extracting the desired top n values defined by another matrix and then dividing by that top n value to get the average at each element, but also found that there were differences between 1 and 2 as seen in your example.

1. sorted_ivt and sqrt(sortx.^2+sorty.^2) (equal here)

2. cummean_sorted_ivt and sqrt(cummeanx.^2+cummeany.^2)

[Walter Roberson]

"Is there any way to reduce this difference?"

You are doing the equivalent of assuming that

A*f(x1, y1) + B*f(x2, y2) + C    %linear combination of values

should equal f(A*x1+B*x2 + C, A*y1+B*y2 + C) %linear combination of positions

for f(x,y) = sqrt(x^2+y^2)

If we simplify with C = 0, and B = 1 - A (that is, that A is the proportion of the way between the two locations), then with that function we can demonstrate that the maximum difference is at A = 1/2, and then the critical points are at A = 1x2 = -x1 and y2 = -y1, which yields a difference between the two functions that works out to sqrt(x1^2 + y1^2) . However, this difference can be made arbitrarily large just by moving x1 and y1. Note that the interpolated position being evaluated at would work out to be f(0,0) -- which is to say that if you allow me to pick the second point as being the reflection of the first point in the X and Y axes, then the point of evaluation would be (0, 0) and if the function is 0 at (0, 0) and your function is the same no matter what the sign of x and y are, such as for your function, then clearly the difference in values would work out to be the same as the function value itself.

For any other point, for that IVT, assuming evaluation half way along the line joining the points, then the difference is

(1/2)*(x1^2+y1^2)^(1/2) + (1/2)*(x2^2+y2^2)^(1/2) - (1/2)*(x1^2+2*x1*x2+x2^2+(y1+y2)^2)^(1/2)

which is definitely not 0.

How can you reduce the difference? By not allowing the second point to be that far from the first point. Which will probably require that you give up your current method of selecting points.

[Vesp]

Thank you, this has been so insightful. I will have to think of another method of selecting the points. I appreciate all of your patience, time and help.