Random Integers from a preconstructed non uniform distribution

Hi, I am looking to select random integers from a specified range that is not uniform. randi([1,interval]) does the job but this is a uniform distribution whereas I wish to change the probability mass function of the distribution on each iteration, the pmf is calculated with accumarray. I apologise if this is unclear but have posted the exact code below for clarity.

x=[1,2];
for j=1:T
    if p/(p+j)>rand

        x=[x,x(end)+1];

    else  

        pmf = accumarray(x(2:length(x))',1); pmf = pmf./sum(pmf);
        x=[x,x(randi([2,x(end)]))];

Many thanks

[Walter Roberson]

It appears to me that you are changing the pmf according to what has already been generated ??

[Michael]

Yes, the pmf and therefore probability of selecting a certain value from 'x' should change on every iteration, e.g. the more 2's there are in 'x' the higher chance we have of selecting 2 to be added into 'x'. Thanks for responding.

Answer 1

[Muthu Annamalai]

A well known approach is called the 'Box Mueller' method. You generate a uniform random variable [0,1] to index into the cumulative distribution function (of the target PDF), and then use it to lookup the corresponding probability interval.

You are access the PDF by inverting the CDF.

As usual Wikipedia explains in detail, http://en.wikipedia.org/wiki/Inverse_transform_sampling_method.

[Al in St. Louis]

That doesn't answer the question since there is no PDF. There is a PMF. [I was wrong. See John D'Errico's comment below.]

[John D'Errico]

Actually, it does. A PMF is essentially the same thing, except it is discrete. You just need to know how to do the lookup. You should notice that the word INDEX was used in the answer.

The simplest way to do this in current MATLAB is to use the discretize function, however, I could have done it using several other tools with not much effort.

X = [2 3 5 7];
P = [.1 .2 .4 .3]; % probabilities for each element in X
N = 100;
Y = X(discretize(rand(1,N),cumsum([0,P])));

[Al in St. Louis]

Thanks! I wasn't getting it when I commented. I was here because I needed to create a nonuniform distribution over the nonnegative integers. At some point later, I realized that I could build a CDF for a PMF. I ended up using interp1([0 cdf.y],[-1 cdf.x],rand(1,N),'next'). I guess it's not elegant, but it worked for me. I am not familiar with the discretize function. The cdf.x vector starts with zero and goes to an integer large enough that the CDF is within 1e-14 of 1. The cdf.y vector starts with the probability of 0 and has the cumulative sum of the PMF for each value in cdf.x.

[John D'Errico]

interp1 was actually one of the alternative methods I could have offered. I even considered suggesting it as such when i said there were other ways. The only flaw with using interp1 is it is perhaps less obvious why that solution works, and when it is necessary to debug code next month or next year, someone might find themselves scratching their head over why you used an interpolation tool there. I think discretize may be slightly more clear to a reader of the code.

Why? In my eyes, interp1 pushes the user to think of the problem in a continuous sense, which is what interpolation is classically used to solve. (Then the careful reader will see that interp1 was used with a discrete method at the end.) However, discretize pushes the user to think of the problem in a fully discrete sense. As such, it takes less of a mental leap to see what was done.

Regardless, I ALWAYS strenuously recommend documenting all code blocks. A line of comment is free, and if that is all that was needed to explain what, how, and why something was done as it was, that can save some brain sweat later on.

Anything that works is good of course.