Thread Subject:

generate random numbers based on probability of occurence

Hi Masters~

I am trying to solve a problem like this: my simulation need a wind from certain direction with certain speed (so, two numbers needed at one time). it must base on a probability of occurrence of the wind.
There are 36 different directions, 10 degree each direction.

for direction from 0 to 260, for each 10 degree the Probability(8m/s) = 0.004
                                                                        the Probability(12m/s) = 0.008
                                                                        the Probability(17m/s) = 0.011
for direction of 270 and 360, the Probability(8m/s) = 0.004
                                         the Probability(12m/s) = 0.010
                                         the Probability(17m/s) = 0.014
for direction of 280 and 350, the Probability(8m/s) = 0.004
                                         the Probability(12m/s) = 0.012
                                         the Probability(17m/s) = 0.013
...and for other directions too.

until now i have no idea how to generate this wind, both direction and speed...

"NUS" wrote in message <jkvek1$5b0$1@newscl01ah.mathworks.com>...
> Hi Masters~
>
> I am trying to solve a problem like this: my simulation need a wind from certain direction with certain speed (so, two numbers needed at one time). it must base on a probability of occurrence of the wind.
> There are 36 different directions, 10 degree each direction.
>
> for direction from 0 to 260, for each 10 degree the Probability(8m/s) = 0.004
> the Probability(12m/s) = 0.008
> the Probability(17m/s) = 0.011
> for direction of 270 and 360, the Probability(8m/s) = 0.004
> the Probability(12m/s) = 0.010
> the Probability(17m/s) = 0.014
> for direction of 280 and 350, the Probability(8m/s) = 0.004
> the Probability(12m/s) = 0.012
> the Probability(17m/s) = 0.013
> ...and for other directions too.
>
> until now i have no idea how to generate this wind, both direction and speed...
- - - - - - - - - - -
  Use matlab's 'rand' function. Prepare a 36 by 3 array of probabilities in which each of the 36 rows corresponds to a direction, each of the 3 columns corresponds to a speed, and each entry gives the probability of the corresponding combination. The sum of all 108 entries should be an exact 1.000 . Let P be such a 36-by-3 array. Then do this:

 P = P/sum(P(:)); % Just to make sure their sum is 1.000
 C = [0;cumsum(P(:))]; C(end) = 2;
 R = rand(n,1); % Generate n random numbers between 0 and 1
 [~,b] = histc(R,C);
 D = mod(b-1,36)+1;
 S = ceil(b/36);

Then D and S will be column vectors of n elements each. The elements of D are indices from 1 to 36 corresponding to the 36 directions. The elements of S are indices from 1 to 3 corresponding the the 3 possible speeds. Each combination will occur with the probability that you placed in its slot in P.

Roger Stafford

Thank you Roger so much for your reply.

i tried, but still got some errors. And i am not really understand how the random generation uses the probability.

"Roger Stafford" wrote in message <jkvo9j$a2l$1@newscl01ah.mathworks.com>...
> "NUS" wrote in message <jkvek1$5b0$1@newscl01ah.mathworks.com>...
> > Hi Masters~
> >
> > I am trying to solve a problem like this: my simulation need a wind from certain direction with certain speed (so, two numbers needed at one time). it must base on a probability of occurrence of the wind.
> > There are 36 different directions, 10 degree each direction.
> >
> > for direction from 0 to 260, for each 10 degree the Probability(8m/s) = 0.004
> > the Probability(12m/s) = 0.008
> > the Probability(17m/s) = 0.011
> > for direction of 270 and 360, the Probability(8m/s) = 0.004
> > the Probability(12m/s) = 0.010
> > the Probability(17m/s) = 0.014
> > for direction of 280 and 350, the Probability(8m/s) = 0.004
> > the Probability(12m/s) = 0.012
> > the Probability(17m/s) = 0.013
> > ...and for other directions too.
> >
> > until now i have no idea how to generate this wind, both direction and speed...
> - - - - - - - - - - -
> Use matlab's 'rand' function. Prepare a 36 by 3 array of probabilities in which each of the 36 rows corresponds to a direction, each of the 3 columns corresponds to a speed, and each entry gives the probability of the corresponding combination. The sum of all 108 entries should be an exact 1.000 . Let P be such a 36-by-3 array. Then do this:
>
> P = P/sum(P(:)); % Just to make sure their sum is 1.000
> C = [0;cumsum(P(:))]; C(end) = 2;
> R = rand(n,1); % Generate n random numbers between 0 and 1
> [~,b] = histc(R,C);
> D = mod(b-1,36)+1;
> S = ceil(b/36);
>
> Then D and S will be column vectors of n elements each. The elements of D are indices from 1 to 36 corresponding to the 36 directions. The elements of S are indices from 1 to 3 corresponding the the 3 possible speeds. Each combination will occur with the probability that you placed in its slot in P.
>
> Roger Stafford

"NUS" wrote in message <jl0qcj$ls6$1@newscl01ah.mathworks.com>...
> i tried, but still got some errors. And i am not really understand how the random generation uses the probability.
- - - - - - - -
  Can you please show me the errors you received?

Roger Stafford

I am quite new on matlab programming.. but there are other cases with defined direction and wind speed in the program which works fine. with the part to generate random direction and speed, here are some problems i dont know why...

my code for this part is like this:
----------------------------------------------
%%generate wind direction and speed
        P = zeros(36,3,'double'); % matrix of 36 by 3, means 36 directions,
                                                % with 3 possible wind speed in each direction
        for a = 1:1:27
            P(a,1) = 0.004;
            P(a,2) = 0.009;
            P(a,3) = 0.011;
        end
        P(28,1) = 0.004; P(28,2) = 0.01; P(28,3) = 0.0135;
        P(29,1) = 0.004; P(29,2) = 0.012; P(29,3) = 0.0165;
        P(30,1) = 0.004; P(30,2) = 0.014; P(30,3) = 0.0195;
        P(31,1) = 0.004; P(31,2) = 0.013; P(31,3) = 0.0315;
        P(32,1) = 0.004; P(32,2) = 0.019; P(32,3) = 0.037;
        P(33,1) = 0.004; P(33,2) = 0.013; P(33,3) = 0.0315;
        P(34,1) = 0.004; P(34,2) = 0.014; P(34,3) = 0.0195;
        P(35,1) = 0.004; P(35,2) = 0.012; P(35,3) = 0.0165;
        P(36,1) = 0.004; P(36,2) = 0.01; P(36,3) = 0.0135;
            
        P = P/sum(P(:)); % Just to make sure their sum is 1.000
        C = [0;cumsum(P(:))];
        C(end) = 2;
        R = rand(n,1); % Generate n random numbers between 0 and 1
        [~,b] = histc(R,C);
        D = mod(b-1,36)+1; %direction
        S = ceil(b/36); %speed choice
        
        theta = D*10 - 10;
        speed = zero(1,3); speed(1) = 8; speed(2) = 12; speed(3) = 17;
        U0 = speed(S);

----------------------------------------------
the error is: ??? Undefined function or variable 'n'.

Error in ==> uwflo at 181
        R = rand(n,1); % Generate n random numbers between 0 and 1

Error in ==> validate>@(x)fitness(x,FitnessFcnArgs{:}) at 136
    fitness = @(x) fitness(x,FitnessFcnArgs{:});

Error in ==> fcnvectorizer at 14
            y(i,:) = feval(fun,(pop(i,:)));

Error in ==> makeState at 47
        Score =
        fcnvectorizer(state.Population(initScoreProvided+1:end,:),FitnessFcn,1,options.SerialUserFcn);
        
Error in ==> galincon at 18
state = makeState(GenomeLength,FitnessFcn,Iterate,output.problemtype,options);

Error in ==> ga at 282
        [x,fval,exitFlag,output,population,scores] = galincon(FitnessFcn,nvars, ...

Error in ==> genetic at 14
[x,fval,exitflag,output,population,score] = ...

Error in ==> main at 31
    [C C1 C2 C3 C4] = genetic(M,lb,ub,200,30000,@uwflo);

Error in ==> run at 74
    evalin('caller',[script ';']);

Caused by:
    Failure in user-supplied fitness function evaluation. GA cannot continue.

On 3/29/2012 11:47 AM, NUS wrote:
...

> P = P/sum(P(:)); % Just to make sure their sum is 1.000
> C = [0;cumsum(P(:))]; C(end) = 2;
> R = rand(n,1); % Generate n random numbers between 0 and 1
> [~,b] = histc(R,C);
> D = mod(b-1,36)+1; %direction
> S = ceil(b/36); %speed choice
> theta = D*10 - 10;
> speed = zero(1,3); speed(1) = 8; speed(2) = 12; speed(3) = 17;
> U0 = speed(S);
>
> ----------------------------------------------
> the error is: ??? Undefined function or variable 'n'.
...

Well even in your comment you say you're going to generate 'n' rngs but
you didn't define 'n' just like the error says...

--

-.- I put 2 to replace n, Now the code is working perfectly.

Thank you so much Roger :)

it is not working again...i dont know why.

    Yi = coord(:,1);
    Xi = coord(:,2);
......
       P = P/sum(P(:)); % Just to make sure their sum is 1.000
        C = [0;cumsum(P(:))];
        C(end) = 2;
        R = rand(2,1); % Generate 2 random numbers between 0 and 1
        [~,b] = histc(R,C);
        D = mod(b-1,36)+1; %direction
        S = ceil(b/36); %speed choice
        
        theta = D*10 - 10;
        speed = zeros(1,3,'int8'); speed(1) = 8; speed(2) = 12; speed(3) = 17;
        U0 = speed(S);
                    
        
        %%change coordinates
        x = cos(theta)*Xi - sin(theta)*Yi;
        y = sin(theta)*Xi + cos(theta)*Yi;

=========
now the error is ??? Error using ==> mtimes
Inner matrix dimensions must agree.
Error in ==> uwflo at 192
        x = cos(theta)*Xi - sin(theta)*Yi;
======================
but theta, Xi, Yi are all non-matrix. why this error would occur?


"NUS" wrote in message <jl2ak3$qst$1@newscl01ah.mathworks.com>...
> -.- I put 2 to replace n, Now the code is working perfectly.
>
> Thank you so much Roger :)

On 3/29/2012 2:12 PM, NUS wrote:
> it is not working again...i dont know why.
>
> Yi = coord(:,1);
> Xi = coord(:,2);
> .......
> P = P/sum(P(:)); % Just to make sure their sum is 1.000
> C = [0;cumsum(P(:))]; C(end) = 2;
> R = rand(2,1); % Generate 2 random numbers between 0 and 1
> [~,b] = histc(R,C);
> D = mod(b-1,36)+1; %direction
> S = ceil(b/36); %speed choice
> theta = D*10 - 10;
> speed = zeros(1,3,'int8'); speed(1) = 8; speed(2) = 12; speed(3) = 17;
> U0 = speed(S);
> %%change coordinates
> x = cos(theta)*Xi - sin(theta)*Yi;
> y = sin(theta)*Xi + cos(theta)*Yi;
> =========
> now the error is ??? Error using ==> mtimes
> Inner matrix dimensions must agree.
> Error in ==> uwflo at 192
> x = cos(theta)*Xi - sin(theta)*Yi;
> ======================
> but theta, Xi, Yi are all non-matrix. why this error would occur?

No, they're not; according to the code above they're two vectors--the
first and second columns of coord. So, if there are more than two
values in coord then Xi and Yi will be vectors of length size(coord,1)

Now if you meant for Xi, Yi to be a single entry from coord, you would
need the reference in a loop and Xi=coord(idx,1) where idx is the
desired row in the coord matrix/array and similar for Yi

--

"NUS" wrote in message <jl2c6b$2vp$1@newscl01ah.mathworks.com>...
> it is not working again...i dont know why.
> Yi = coord(:,1);
> Xi = coord(:,2);
- - - - - - - -
  In answer to your question, "i am not really understand how the random generation uses the probability", the C = [0,cumsum(P(:))] vector sets up 180 different subintervals from 0 to 1, in which each interval's length is equal to the corresponding probability in P. This means that the probability that a 'rand' number will lie in each particular interval will be just that value in P. The 'histc' function detects which intervals the 'rand' values have fallen in and places their corresponding linear indices, 1 to 108, in 'b' (bin number.) The 'mod' and 'ceil' operations convert from these linear indices back to the corresponding subscripts in P, which will be the direction and speed indices. I hope this abbreviated explanation is clear.

  I neglected to tell you that I used 'n' to indicate how many random numbers you want to generate. It is most efficient to generate them all in one call to 'rand' as I did in this code and save the corresponding wind directions and speeds in 'theta' and 'speed' arrays to use in your simulation afterwards.

  However that means that your 'theta' and 'speed' quantities will be vectors, not scalars, and must be treated accordingly. Since your Xi and Yi quantities are also vectors, you must decide how to handle them properly. Your most recent error message is telling there is something wrong with the way you were doing this. I can't advise you on that aspect of your simulation. I don't know how the various different coordinates are combined with random wind directions and speeds. You either have to use some kind of for-loop method or perhaps some vectorized method.

  One final note of caution. Unless the number 'n' is an exceedingly large number, you must not expect that the actual observed distribution among your 180 possibilities will be just the same as the given probabilities. These probabilities will only be the mean statistical values over a large number of observations. Quite often in this newsgroup users have made comparatively small numbers of Monte Carlo trials and then complained that their percentages were not in accord with the required probabilities. However stochastic processes like yours do not act like that. They require very large numbers of observations to begin to approach their statistical mean values. For example, there is less than one chance in four that a penny flipped ten times will produce exactly five heads and five tails, even though those are the expected numbers. With 180 different possibilities you will need
very, very much larger numbers than ten.

Roger Stafford

Thank you, it is so clearly written.
so, according to your explanation, then the 2-demension array for P should be 3*36, not 36*3, right?
if i use the 36*3 array, i need to change the code to:

direction = ceil(b/3);
speed = mod(b-1, 3)+1;

am I right?

what is the meaning of C(end) = 2 ?why the histc(x,edgexs), edges is C? why....what is [~,b] = histc(R,C)? does all these return n randome values to b vector?


"Roger Stafford" wrote in message <jl2jin$ro3$1@newscl01ah.mathworks.com>...
> "NUS" wrote in message <jl2c6b$2vp$1@newscl01ah.mathworks.com>...
> > it is not working again...i dont know why.
> > Yi = coord(:,1);
> > Xi = coord(:,2);
> - - - - - - - -
> In answer to your question, "i am not really understand how the random generation uses the probability", the C = [0,cumsum(P(:))] vector sets up 180 different subintervals from 0 to 1, in which each interval's length is equal to the corresponding probability in P. This means that the probability that a 'rand' number will lie in each particular interval will be just that value in P. The 'histc' function detects which intervals the 'rand' values have fallen in and places their corresponding linear indices, 1 to 108, in 'b' (bin number.) The 'mod' and 'ceil' operations convert from these linear indices back to the corresponding subscripts in P, which will be the direction and speed indices. I hope this abbreviated explanation is clear.
>
> I neglected to tell you that I used 'n' to indicate how many random numbers you want to generate. It is most efficient to generate them all in one call to 'rand' as I did in this code and save the corresponding wind directions and speeds in 'theta' and 'speed' arrays to use in your simulation afterwards.
>
> However that means that your 'theta' and 'speed' quantities will be vectors, not scalars, and must be treated accordingly. Since your Xi and Yi quantities are also vectors, you must decide how to handle them properly. Your most recent error message is telling there is something wrong with the way you were doing this. I can't advise you on that aspect of your simulation. I don't know how the various different coordinates are combined with random wind directions and speeds. You either have to use some kind of for-loop method or perhaps some vectorized method.
>
> One final note of caution. Unless the number 'n' is an exceedingly large number, you must not expect that the actual observed distribution among your 180 possibilities will be just the same as the given probabilities. These probabilities will only be the mean statistical values over a large number of observations. Quite often in this newsgroup users have made comparatively small numbers of Monte Carlo trials and then complained that their percentages were not in accord with the required probabilities. However stochastic processes like yours do not act like that. They require very large numbers of observations to begin to approach their statistical mean values. For example, there is less than one chance in four that a penny flipped ten times will produce exactly five heads and five tails, even though those are the expected numbers. With 180 different possibilities you will need

> very, very much larger numbers than ten.
>
> Roger Stafford

"NUS" wrote in message <jl37qh$634$1@newscl01ah.mathworks.com>...
> Thank you, it is so clearly written.
> so, according to your explanation, then the 2-demension array for P should be 3*36, not 36*3, right?
> if i use the 36*3 array, i need to change the code to:
>
> direction = ceil(b/3);
> speed = mod(b-1, 3)+1;
>
> am I right?
>
> what is the meaning of C(end) = 2 ?why the histc(x,edgexs), edges is C? why....what is [~,b] = histc(R,C)? does all these return n randome values to b vector?
- - - - - - - - -
  I will attempt to further explain the use of 'histc' here. To simplify matters, suppose P is a vector with only four elements: P = [0.29;0.13;0.37;0.21]. Then C would be C = [0;0.29;0.42;0.79;2.00]. If we use this as the 'edges' vector in 'histc', the probability of a 'rand' value falling between 0 and 0.29 is just 0.29-0=0.29, between 0.29 and 0.42 it is 0.42-0.29=0.13, between 0.42 and 0.79 it is 0.79-0.42=0.37, and finally between 0.79 and 2.00 it is 1.00-0.79=0.21 (remember that rand must always be less than 1.00) These are precisely the required probability values in P. The second argument returned by 'histc' gives the corresponding bin numbers for these four intervals and therefore amounts to the linear indices into P for the 'rand' numbers.

  The reason I changed the last element in C to 2 is that there is a very tiny chance that due to round off errors in computing the cumsum(P(:)) its final number might be smaller than 1.00 and a rand value might be generated greater than or equal to it (actually this is a most unlikely happening.) However if it did happen, you would get an index of 109 or 0 which would mess up things. For some incomprehensible reason I chose to put a 2 there. I could just as well have set it to 1 or +inf. If you like, you can put a 1 there instead which is very close to what is already there.

  As to changing to 3-by-36, that all depends on what you are going to do with Xi and Yi. There wasn't anything in my explanation itself that would require a 3-by-36 change. Linear indexing would work equally well with either 36-by-3 or 3-by-36. In case it is necessary to make the change, your modifications of 'direction' and 'speed' look correct to me.

Roger Stafford

Try these:

http://en.wikipedia.org/wiki/Inverse_transform_sampling
http://www.mathworks.com/matlabcentral/fileexchange/7309-randraw

hi there,
i have a question for Roger Stafford. i'm doing something quite similar but with wave data instead. My 'P' is a 19x21 matrix containing occurence probabilites of certain wave conditions. where rows represent wave heights (D) and columns represent wave periods (S). let's call them D and S respectively so as to avoid confusion in the current thread. My question is about how you obtain the indices.

currently this section of my code looks like this:

 D = mod(b-1,19)+1;
 S = ceil(b/19)

why do you have '+1' to D?
should i change the '19' in S to '21' or should i keep it as it is?

sorry, i don't quite understand how these two lines work as i'm fairly new to matlab as well.

thanks

"Joshua " <s0567996@ed-alumni.net> wrote in message <kh0ik4$n6v$1@newscl01ah.mathworks.com>...

I'm not Roger. But I will asnwer

> why do you have '+1' to D?

to make D varying from 1 to 19

> should i change the '19' in S to '21' or should i keep it as it is?
>
No 19 is OK

Alternatively you can compute D and S from b by

[D S] = ind2sub(size(P),b)

Bruno

hi bruno, thanks very much!

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <kh1kp8$lb0$1@newscl01ah.mathworks.com>...
> "Joshua " <s0567996@ed-alumni.net> wrote in message <kh0ik4$n6v$1@newscl01ah.mathworks.com>...
>
> I'm not Roger. But I will asnwer
>
> > why do you have '+1' to D?
>
> to make D varying from 1 to 19
>
> > should i change the '19' in S to '21' or should i keep it as it is?
> >
> No 19 is OK
>
> Alternatively you can compute D and S from b by
>
> [D S] = ind2sub(size(P),b)
>
> Bruno