Mean of multiple large matrices with NaN's

I have 23 individual matrices each are 385x781x365. They are saved as .mat files, each are ~100mb. I would like to take the mean of all 23 along the third dimension. Unfortunately, there is too much data to load all 23, concatenate them, and use mean(concatenated_array,3).

As an alternative, I have resorted to a loop (dummy data for your use):

 data = randn(385,781,365,23); 
 for n = 1:23
    if n == 1
        data_mean = data(:,:,:,n)./23;
    else
        data_mean = data_mean + data(:,:,:,n)./23;
    end
 end

 data_mean2  =   squeeze(mean(data,4));

If you look at the data, you will see that data_mean and data_mean2 are equivalent (isequal will return 0 because of rounding, though).

So, this works ok. However, a problem arises because the 23rd matrix is actually 385x781x274. Thus, on the final iteration the loop fails because the matrices are not the same size.

I have tried setting the 275:365 of the final matrix equal to the mean:

 data = randn(385,781,365,23); 
 for n = 1:23
    if n == 1
        data_mean = data(:,:,:,n)./23;
    elseif n == 23
        data(:,:,275:365,n)     =   data_mean(:,:,275:365);
        data_mean               =   data_mean + data(:,:,:,n)./23;
    else
        data_mean = data_mean + data(:,:,:,n)./23;
    end
 end

However you will notice that this is not right.

Ideally, I would like to set the 275:365 of the final matrix to nan, but then upon division the entire data_mean becomes NaN for 275:365.

Any ideas how to fix? Ideally a solution that does not require the loop. Thank you.