0

i'm developing a program in octave that i will explain as i put the code. So i have this matrix in a file called matprec.m:

function [res1] = matprec()
   res1 = [
1,2001,1,2,0.00;
1,2001,1,5,5.33;
2,2001,1,5,4.57;
3,2001,1,5,5.33;
4,2001,1,5,5.59;
5,2001,1,5,4.32;
2,2001,1,13,0.00;
3,2001,1,13,0.00;
4,2001,1,13,0.00;
3,2001,1,30,30.73;
2,2001,2,1,1.02;
3,2001,2,1,1.52;
4,2001,2,1,1.78;
5,2001,2,1,1.27;
1,2001,2,2,1.78;
2,2001,2,2,1.27;
3,2001,2,2,1.78;
4,2001,2,2,2.03;
5,2001,2,2,1.78;
1,2001,3,4,18.03;
3,2001,3,4,15.75;
5,2001,3,4,17.53;
1,2001,3,5,13.46;
2,2001,3,5,12.19;
3,2001,3,5,11.94;
4,2001,3,5,9.65;
5,2001,3,5,10.92;
2,2001,4,30,0.00;
4,2001,4,30,0.00];
format short g
return
endfunction

So in this matrix the first column is just the station where we measure the amount of precipitation, the second is the year, the third is the month, the fourth is the day and the fifth is the value of precipitation. And what i want to do in another file is call this matrix and do the following calculus, in the month 1 i want do the average on all the days for example:

in month 1 day 5 i have 5 values 5.33, 4.57, 5.33, 5.59, 4.32, so i would do

(5.33 + 4.57 + 5.33 + 5.59 + 4.32)/5 = 5.028

And i want to do that for all the days and when i have all the days i would add them all to know the amount of precipitation in that month, and do that for all the 4 months.

That's how the program should function and i already accomplished that with this piece of code:

Result        = matprec();
month1Indices = Result(:,3) == 1;
month1Rows    = Result(month1Indices, :);
day5Indices   = month1Rows(:,4) == 5;
day5Rows      = month1Rows(day5Indices , :);
  mean(day5Rows(:,5));

Result2        = matprec();
month1Indices2 = Result2(:,3) == 1;
month1Rows2    = Result2(month1Indices2, :);
day5Indices2  = month1Rows2(:,4) == 30;
day5Rows2      = month1Rows2(day5Indices2 , :);
  mean(day5Rows2(:,5));

  mes1 = mean(day5Rows(:,5)) + mean(day5Rows2(:,5));

lol = [mes1, mes2, mes3, 0, 0, 0, 0, 0, 0, 0, 0, 0];
  x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12];
  plot(x,lol);

But that only work for the first month, so i'm trying to do a while loop to do that for all the month's and all the day's and plot it, i have this but it just continue forever:

a = 1
b = 1
while (a <= 4 && b <= 30)
Result        = matprec();
month1Indices = Result(:,3) == a;
month1Rows    = Result(month1Indices, :);
day5Indices   = month1Rows(:,4) == b;
day5Rows      = month1Rows(day5Indices , :);
  mean(day5Rows(:,5))
  x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12];
  plot(x,mean(day5Rows(:,5)));
  endwhile

Sory i know it's a bit long, but i need to explain it, thanks.

2 Answers2

0

I rewrote this a bit, since i think you overcomplicating things by using so many different variables and hardcoding of numeric values. Often it's easier to extract just exactly what you need:

months = res1(:, 3);
prec = res1(:, 5);

m_start = 1; # month where to start averaging
m_end = 4;   # month where to end averaging

# prepare output vector
avg = zeros(m_end - m_start + 1, 1);

for m = m_start:m_end
    temp = prec(months == m);
    avg(m) = mean(temp);
end

plot(avg)
grid

In your example, you only used the non-zero values, is this intended? If you want this behaviour, i suggest using a small epsilon value as threshold:

for m = m_start:m_end
    temp = prec(months == m);
    # --------------------------
    # ignore smaller values
    eps = 1e-12;
    temp = temp(temp > eps);
    # --------------------------
    avg(m) = mean(temp);
end
pschulz
  • 1,406
  • 13
  • 18
0

The calculation you describe is a little confusing. I'll assume that the first step is to determine the mean precipitation on each day in a month (averaged over the stations)

mean_precipitation = accumarray([res1(:,3), res1(:,4)], res1(:,5), [], @mean);

Each row of the matrix corresponds to a month, and each column is a day. The accumarray function does the majority of the work that you are doing in a loop. The first argument [res1(:,3), res1(:4)] is a two column matrix, where the first column specifies the month (row) and the second column specifies the day (column) that each 'accumulated' result belongs to. The second argument res1(:,5) are the values which are 'accumulated' into each final result, and the final argument @mean is the type of accumulating function to apply to the collection of each month-day's data. See help accumarray for an alternative description (or google) - it can be a bit confusing how it works at first.

There are two alternative interpretations of the next step in your calculation.

Mean precipitation on a rainy day in a month

If you want the average precipitation of days where there was precipitation over a month then:

sum(mean_precipitation, 2) ./ sum(mean_precipitation > 0, 2)

ans =

   17.8790
    1.5627
   14.3677
       NaN

This will give a NaN for the final month as there were no days with precipitation recorded, so the mean is not defined.

Mean precipitation on all days in a month

If you want the average precipitation over all days in the month, then you need to work out the number of days in each month. This can be found via eomday() (I assume that all observations are in year 2001), the calculation is then

sum(mean_precipitation, 2) ./ eomday(2001, (1:size(mean_precipitation, 1))')

ans =

   1.15348
   0.11162
   0.92695
   0.00000

Hopefully this gives you some ideas about how to perform your calculation.

Why no while loops?

While loops are an older imperative style of writing code. Sometimes they are useful for clarity, but sometimes they are slower than using established functions which may 'vectorise' the calculation (see Why does vectorized code run faster than for loops in MATLAB?)

Community
  • 1
  • 1
stephematician
  • 844
  • 6
  • 17