Firstly, I don't have a stats background, so please accept my apologies for any errors or misunderstandings in the question below.

I'm trying to use R to draw values from a discrete probability distribution across a defined range of integers, e.g. 0:100. In my application these integers represent balls, labelled 0 to 100, which are to be drawn with replacement from a bag.

I'm looking for a discrete distribution that I can draw from which I am able to vary (perhaps by modifying one or more parameters) such that it starts off equivalently to a discrete uniform distribution (with every ball having the same probability of being drawn), and becomes increasingly left-skewed as a parameter is modified (i.e. balls with higher numbers have a greater and greater probability of being drawn).

The idea is that in its initial state, the distribution would be completely uniform, and as a particular parameter is adjusted it becomes increasingly left-skewed. Or, in practical terms, the bag would start with the same amount of balls of each number, and on each round there would be more and more balls of higher numbers, and fewer of lower numbers.

i.e. starting with this:


and going to something like this:

Left skewed

Is there a discrete distribution that has these properties, and if so, how can I use its probability mass function in R to determine the probability of each ball being drawn for any one such distribution?


2 Answers 2


You could look at the beta-binomial distribution, $X\sim\text{BB}(n,\alpha,\beta)$. If $\alpha=\beta=1$, then it is a discrete uniform distribution, if $\alpha\ge 1$ and $\beta<1$ then it is a discrete left-skewed distribution.

In R you could use the extraDistr package. For example:

> library(extraDistr)
> x <- 0:100
> y <- dbbinom(x, 100, 1, 1)
> plot(x, y, type="h")

enter image description here

> y <- dbbinom(x, 100, 1, 0.9)
> plot(x, y, type="h")

enter image description here

> y <- dbbinom(x, 100, 1.1, 0.9)
> plot(x, y, type="h")

enter image description here

You can draw (pseudo)random numbers from a beta-binomial distribution by rbbinom.

  • $\begingroup$ Beautiful, this is exactly what I needed. Thank you so much! $\endgroup$ Aug 10, 2020 at 7:06

Many simulation problems seek a large random sample of independent and identically distributed (IID) random observations, all from the same population. The built-in simulation procedures in R are mainly designed with such IID samples in mind.

So it is necessary to do some extra work to change the distribution during the sampling process.

100 Coin Tosses: Suppose you're tossing a fair coin 100 times independently. One hundred Bernoulli trials with $P(H) = 1/2.$

You could use the following for 100 IID tosses all with P(H) = 0.5 and then plotting results in sequence. (The last argument of plot causes individual data points to be connected by lines.)

x = rbinom(100, 1, .5)
plot(x, ylim=c(-.3,1.3), pch=20, type="b")

enter image description here

Heads (at 1) are about equally frequent at beginning and end of my 100 tosses.

Now suppose your coin somehow becomes increasingly more biased in favor of Heads as the process goes through the successive tosses. You could get this result by specifying different Heads probabilities for each of the 100 tosses.

p = seq(.5, 1, len = 100)
y = rbinom(100, 1, p)
plot(y, ylim=c(-.3,1.3), pch=20, type="b")

enter image description here

Drawing 100 chips from urn without replacement. Suppose I have an urn with 100 chips numbered from 1 through 100. I will draw them out one at a time without replacement until they are all gone. Even if I draw chips at random, this is not an independent process, each time I remove a chip, that affects my choices on the next draw, and finally on the last draw I can only take the one remaining chip. But by drawing the chips so that each of the remaining chips is equally likely to be chosen at each step, I have essentially put drawn chips into a random order.

In the sample function the first argument is the (initial) population and the second is the number of draws.

w = sample(1:100, 100)
plot(w, ylim=c(-.3,100.3), pch=20, type="b")

enter image description here

At the beginning, and still at the end, of the 100 draws I have gotten a mixture of low and high-numbered chips.

If I want to bias the process to that I tend to get larger numbers first, then I can use a probability vector of length 100 to say at each step that higher-numbered chips have priority. Priority numbers for missing chips are ignored for each draw and R re-scales remaining priority numbers to sum to $1$ again after each draw.

priority = 1:100
v = sample(1:100, 100, p=priority)
plot(v, ylim=c(-.3,100.3), pch=20, type="b")

enter image description here

Now it is clear the higher-numbered balls are chosen early in the process, with lower numbers remaining toward the end.

Notes: (1) If you have some more intricate kinds of biasing tendencies in mind, then you could us a hand-written loop (one step for each draw) and change preferences as the process continues. For example, you could increase the span of possible values with each draw from the urn:

m = 100;  u = numeric(m)
for(i in 1:m) {
 u[i] = sample(1:i, 1)  # note population changes
plot(u, ylim=c(-.3,100.3), pch=20, type="b")

enter image description here (2) I have put a set.seed statement at the start of each bloc of code. If you retain them, you will get exactly the same simulations I have shown. If you remove them, then R will pick an unpredictable seed and you will get a fresh simulation demo on each run.

  • $\begingroup$ Hi BruceET, Thanks for your answer. It's really helped me get my head around some of the ways R deals with simulating observations. However, I think I might've been a little unclear in my question: in my specific case the number of balls in the bag isn't necessarily 100, i.e. there could be multiples of each ball number, and I'd be sampling with replacement. Also, I wasn't suggesting that the distribution change after each draw, but rather after each run – so I'd get, e.g., 500 observations from a uniform distribution, then 500 more from a slightly biased one, and so on... $\endgroup$ Aug 10, 2020 at 5:16
  • $\begingroup$ I tried to focus on what would work with built-in R procedures and without writing loops. Of course, if you change the population at each step of a run (as in my Note) or with each run you can do almost anything. You could also program a Markov Chain as in this Q&A. $\endgroup$
    – BruceET
    Aug 10, 2020 at 5:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.