How to pick random numbers from an interval non-uniformly?

I want to pick some numbers (N) from the interval [0,2Pi] randomly. I want all the real number in the [pi/2,pi] have probability P1 and all the others ([0,pi/2]+[pi/2pi]) have probability P2<P1. basically, the sample space has two portions: the individuals within each portion have equal probability (uniform distribution for each portion).

can I extend the "rand" function somehow to do this?

In addition, these random numbers represent "angels" of the points that should be placed on the perimeter of a circle with radius R. IF I assume uniform distribution for the whole sample space [0,2pi], is the following a correct way to choose N sample points:

    angels=rand(N,1);
    angels=0+(2*pi-0)*angels;

'cause most of the times the points seem to be more densely placed in some regions.

Any help would be appreciated

tnx in advance

Answer 1

[Image Analyst]

Basically you want to do inverse transform sampling, as explained in http://en.wikipedia.org/wiki/Inverse_transform_sampling It involves taking the cumulative distribution function. I've attached a demo for picking numbers from a Rayleigh distribution.

For more distributions you can see RANDRAW: http://www.mathworks.com/matlabcentral/fileexchange/7309-randraw

[ramin bba]

Thanks image Analyst. It is the exact thing I was looking for.

P.S.: I assume you deleted the irrelevant comment. tnx for that as well.

regards,

[Image Analyst]

You're welcome - no problem. Thanks for "Accepting" the answer.

Answer 2

[James Tursa]

One method (assuming P1+P2=1):

A = rand(N,1);
I = A < P1; % indexes of the P1 set
A(I) = A(I)*(pi/(2*P1)); % scale the P1 set to [0,pi/2]
A(~I)=(A(~I)-P1)*(3*pi/(2*(1-P1)))+pi/2; % scale the P2 set to [pi/2,2pi]
A = mod(A+pi/2,2*pi); % shift the sets over to the desired final locations

For your circle question, yes you can simply use 2*pi*rand(N,1) to get a uniform distribution of angles in the [0,2pi] range. You can get the same result by using P1=0.25 in the code I supplied. As far as your "densely placed" observation goes, it is likely just due to not many samples being generated. Try it with a large number of samples.

[ramin bba]

tnx a lot, it works! however, I was looking for a more generic form (which is not unfortunately stated in my original question). here is the problem:

choose 5 points in the intervals [pi/4,3pi/4]+[3pi/2,7pi/4].

One difficulty that i am facing is that "5" is not a large number! Is there any way around this?

also, if I want to use your method, do I have to just play around with the last line?

here is how I was thinking of doing it (for the problem I described here):

alpha*[(3pi/4-pi/4)+(7pi/4-3pi/2)]=1 -----> alpha=4/3pi; %probability per unit angel

afterwards, I choose a random number in [0,1]. if it is less than 4/3pi*(3pi/4-pi/4)=2/3; it will be in the first interval. otherwise, it will go into the second one. once I know how many samples I need to choose from each interval, I could simply do a+(b-a)*rand(Nab,1) where a and b are the bounds of a specific interval.

i think this will have the same results as yours, right? (although your code is conciser).

thanks again

[James Tursa]

Can you restate your real question completely? E.g., I don't understand these comments:

choose 5 points in the intervals [pi/4,3pi/4]+[3pi/2,7pi/4].

Choose them how? Uniformly in each? Different probabilities? Or ...?

One difficulty that i am facing is that "5" is not a large number! Is there any way around this?

Why is this a difficulty?

also, if I want to use your method, do I have to just play around with the last line?

The method I gave is specific to the problem you posed. So yes, pretty much everything in the method, including the last line, might have to be adjusted or changed if there is a different problem being solved.

Your last paragraph I don't completely understand, but it seems to imply a for loop of some kind iterated on the samples. This will of course work, but will be slow if you have a lot of samples to process (may not be a concern to you).

[ramin bba]

Thank a lot James for the thorough follow up and sry for the confusion.

here is the actual problem:

Choose N points (N is small, i.e. 5) from some subset(s) of the the interval [0,2pi]?

i.e., assume the subsets are S1 = [pi/4,3pi/4] and S2 = [3pi/2,7pi/4]. Now two case might happen:

- ALL the points in the sample space have equal probability (i.e. a point in S1 is as equally likely to be chosen as a point in S2).

- points in ONE subset have equal probability. for example, the probability of a subset might be proportional to its length (the longer an interval, the more likely it is for its individuals to be chosen. This is why I defined "alpha", probability per unit length, in my previous comment).

regarding the difficulty, as I said the point that I chose seem to be more dense in some regions and you said (in your first reply) that if the number of samples is large, this will not happen! well, the number of sample IS small (N=5). Does this mean that the issue will always happen and there is no way around it?

I hope I was clear this time.

BTW, regarding the totally irrelevant answer below, how can I delete it? I can Only flag it, right?

thanks again for the thorough responses.

regards,

Answer 3

[Sam]

I believe this should do the trick. Just set P2 and P1 to the values you want (their exact value is not important, just their ratio). Essentially what this code does is instead of generating random numbers from 0 to 2*pi, it generates numbers from a larger range (the bigger the P1/P2 ratio, the bigger the range). It then shrinks the region from pi/2 to pi+(P1/P2-1)*pi/2 into the region pi/2 to pi, which increases the density there. Then it just offsets the remainder of the samples greater than pi so that the region from 0 to 2 pi is filled. There might be a way to do this using logical indexing instead of for loops, but I can't think of it now.

P2 = 1;
P1 = 2;
numsamples = 100000;

angles = rand(numsamples,1)*2*pi*(0.75+P1/P2/4);
for count = 1:numsamples,
    if (angles(count) > pi/2) && (angles(count)<(pi+pi*(P1/P2-1)/2))
        angles(count) = pi/2+(angles(count)-pi/2)/(P1/P2);
    elseif angles(count)>(pi+pi*(P1/P2-1)/2)
        angles(count) = angles(count) - pi*(P1/P2-1)/2;
    end
end
hist(angles,200)

[ramin bba]

tnx for the response. A couple of questions and concerns:

- why are you using 0.75 and 1/4 in the 4th line?

- although this would not affect the results, P1+P2 should equal 1 since the two sets comprise the population.

-I am strongly inclined to avoid using "for loops" 'cause I have to do this procedure a lot in my main function.

tnx again

[Sam]

The 0.75 and 1/4 are from the fact that your interval of pi/2 to pi is one quarter of the total [0,2*pi] interval. Similarly, the 3/4 is because the rest is because the rest is 0.75 of the interval. As I said, the absolute values of P1 and P2 are irrelevant so you could set them equal to one with the ratio you want (so you would have P1 = 0.67 and P2 = 0.33 if you want the denser interval to have twice the density).

If you're doing stuff with the distribution, my guess would be that the for loop would not be a major chunk of the computation time. You could do a convoluted logical indexing, but I don't think it would save you much time and it would make the code much less readable. Plus it sounds from you comment below that your number N will be fairly small (like 5). In this case, I don't think you need to worry about computation time for that.