Is there a more efficient way to average multidimensional matrix elements?

I am working with high-dimensional arrays with inconsistent dimensionalities, but all contain 6+ dimensions, and the last 3 always represent 3D space in XYZ. For each x,y,z location, I wish to calculate the mean value for elements that index, collapsing across all other dimensions. I am currently using vectors to gather these values within nested for-loops, and averaging them, as in the following set of code snippets ('betas' is the multidimensional array in question):

First, obtain dimensions of betas:

betasdim=size(betas);

Next, calculate the size of the 3D space, and the number of dimensions I need to collapse over. The last 3 dimensions of betasdim are XYZ:

voxdim=betasdim(length(betasdim)-2:length(betasdim));

Everything else are the dimensions to collapse over (create a vector of colons to collapse over all dimensions other than XYZ:

otherdims=betasdim(1:length(betasdim)-3);
numdims=length(otherdims);
dimwildcard=repmat({':'}, 1, numdims);

Initialize matrix of means.

 meanbetas=repmat([NaN],voxdim);

And finally, the likely inefficient for-loop solution:

for x=1:voxdim(1)
 for y=1:voxdim(2)
  for z=1:voxdim(3)
    voxbetas=betas(dimwildcard{:},x,y,z);%get all beta values for this xyz 
    voxbetas=reshape(voxbetas,1, numel(voxbetas));%reshape to vector
     meanbetas(x,y,z)=nanmean(voxbetas); %average the vector and 
                      %store in corresponding xyz index in the matrix of means
  end
 end
end

Bearing in mind that I need a single value from nanmean() at each index, is there a faster solution other than to iterate through each value of x,y,z?

Answer 1

[Matt J]

 meanbetas=nanmean(reshape(betas,prod(otherdims), []) ,1);

 meanbetas=reshape(meanbetas,voxdim);