Sampling From a Discrete Distribution – How to Perform the Process

distributionsprobabilityrandom-generationweighted-sampling

Assume I have a distribution governing the possible outcome from a single random variable X.
This is something like [0.1, 0.4, 0.2, 0.3] for X being a value of either 1, 2, 3, 4.

Is it possible to sample from this distribution, i.e. generate pseudo random numbers upon each of the possible outcomes given the probability of that outcome. So if I wanted to know what the probability of getting a 2 is, the sampling operation may return 0.34 or something like that.

The reason I ask is that I'm trying to implement an action selection policy for a reinforcement learning method based on a research paper. From what I gather from the paper, the author is able to sample the distribution by "mapping the uniform distribution U[0,1] through cumulative probability density functions obtained by adaptive numerical integration". From this he then samples the transition probabilities for each trial…

I would be grateful for any info on this…

Thanks in advance

Best Answer

Sure. Here's an R function that will sample from that distribution n times, with replacement:

sampleDist = function(n) { 
    sample(x = c(1,2,3,4), n, replace = T, prob = c(0.1, 0.4, 0.2, 0.3)) 
}

# > sampleDist(10)
# [1] 4 2 2 2 2 2 4 1 2 2

If you want to go a little lower level, you can see the actual algorithm used by checking out the R source (written in C):

/* Unequal probability sampling; with-replacement case 
 * n are the lengths of p and perm. p contains probabilities, perm
 * contains the actual outcomes, and ans contains an array of values 
 * that were sampled.
 */

static void ProbSampleReplace(int n, double *p, int *perm, int nans, int *ans)
{
    double rU;
    int i, j;
    int nm1 = n - 1;

    /* record element identities */
    for (i = 0; i < n; i++)
        perm[i] = i + 1;

    /* sort the probabilities into descending order */
    revsort(p, perm, n);

    /* compute cumulative probabilities */
    for (i = 1 ; i < n; i++)
        p[i] += p[i - 1];

    /* compute the sample */
    for (i = 0; i < nans; i++) {
        rU = unif_rand();
        for (j = 0; j < nm1; j++) {
            if (rU <= p[j])
                break;
            }
        ans[i] = perm[j];
    }
}
Related Question