"The data has months 7 to 10 and labelled them as M7, M8, M9, M10 to help me remember the months."
Do not do that! Numbered variables are always a bad idea. Your numbered variables will force you to copy/paste your code for each month just to change the variable name. It will make everything more complicated, not simpler.
"Every M that is recorded for July has a value 7 attributed to it in another row. The same for August with a value 8 in another row etc. I basically indexed for the month == 7 to get the size of that and then extracted the magnitude values for that size in the column with magnitudes"
So you basically have a matrix with at least two columns with one column being the month and the other the magnitude. You may consider using a table for this as tables have plenty of features suited to this type of data. In particular, grouptransform would allow you to do your calculation for all the months at once. Assuming that M is the magnitude column of your matrix and month the month column, your initial extraction of the values corresponding to a month does not work. Assuming that there are N rows corresponding to a month, your code is ignoring which actual rows these are and always extract the first N rows of the matrix. The proper code would be:
thismonthdata = M(month == 7);
I'm a bit confused by your equation "b=loge/(M_ave - M_min)", is by definition 1. I'm assuming your actual equation is , in which case it can be implemented very simply as windowsize = 100
b = log(movmean(M, windowsize, 'EndPoints', 'discard') - movmean(M, windowsize, 'EndPoints', 'discard'));
There's no built-in function to perform your second calculation over a window, so you will have to use a loop. In that case, you may as well do the above in the loop as well, it would go like:
windowsize = 100;
for start = 1:numel(M)-windowsize+1
windowdata = M(start+(1:windowsize));
b = log(mean(windowdata) - min(windowdata));
frequency_table = tabulate(windowdata);
logcumulative_frequency = log10(cumsum(frequency_table(:, 3), 1, 'reverse'));
magnitude = frequency_table(:, 1);
[p, S] = polyfit(magnitude, logcumulative_frequency, 1);
[y_fit, delta] = polyval(p, magnitude, S);
end