How to sample from a discrete distribution? 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 
 A: In python you could do something like
    from scipy.stats import rv_discrete

    x=[1,2,3,4]
    px=[0.1,0.4,0.2,0.3]

    sample=rv_discrete(values=(x,px)).rvs(size=10)

Which would give you 10 samples from the distribution.  You could repeat this then find the proportions of samples that are 2.
A: Yes it is possible and fairly easy, exactly how depends on what tool(s) you are using.
In R it would be sample(1:4, n, prob=c(0.1,0.4,0.2,0.3), replace=TRUE) where n is the number of values you want to sample.
In tools without an equivalent function you can still generate a uniform value and then your RV will equal 1 if it is below 0.1, 2 if it is between 0.1 and 0.5, 3 if between 0.5 and 0.7, and 4 if greater than 0.7 (that is the idea of mapping to the cumulative).
For your example you could also sample uniformly from a set with one 1, four 2's, two 3's, and three 4's to get the same probabilities.
A: 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];
    }
}

A: In Stata: 
In Mata use rdiscrete() as documented at http://www.stata.com/help.cgi?mf_runiform 
In Stata itself, there are various ways. Here's one: 
. gen rnd = runiform()
. gen y = cond(rnd <= 0.1, 1, cond(rnd <= .5, 2, cond(rnd <= .7, 3, 4)))

