Is there a discrete distribution I can use for sampling in R?

Question

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:

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?

Answer 1

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")

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

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

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

Answer 2

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.)

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

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.

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

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.

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

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.

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

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:

set.seed(2020)
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")

(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.