## Using Stata’s random-number generators, part 2, drawing without replacement

Last time I told you that Stata’s **runiform()** function generates rectangularly (uniformly) distributed random numbers over [0, 1), from 0 to nearly 1, and to be precise, over [0, 0.999999999767169356]. And I gave you two formulas,

- To generate
*continuous*random numbers between*a*and*b*, use**generate double u =****(***b***–***a***)*runiform() +***a*The random numbers will not actually be between

*a*and*b*: they will be between*a*and nearly*b*, but the top will be so close to*b*, namely 0.999999999767169356**b*, that it will not matter. - To generate
*integer*random numbers between*a*and*b*, use**generate ui =****floor((***b***–***a***+1)*runiform() +***a***)**

I also mentioned that **runiform()** can solve a variety of problems, including

- shuffling data (putting observations in random order),
- drawing random samples without replacement (there’s a minor detail we’ll have to discuss because
**runiform()**itself produces values drawn*with*replacement), - drawing random samples with replacement (which is easier to do than most people realize),
- drawing stratified random samples (with or without replacement),
- manufacturing fictional data (something teachers, textbook authors, manual writers, and blog writers often need to do).

Today we will cover shuffling and drawing random samples *without* replacement — the first two topics on the list — and we will leave drawing random samples with replacement for next time. I’m going to tell you

- To place observations in random order — to shuffle the observations — type
. generate double u = runiform() . sort u

- To draw without replacement a random sample of
*n*observations from a dataset of*N*observations, type. set seed

*#*. generate double u = runiform() . sort u . keep in 1/*n*I will tell you that there are good statistical reasons for setting the random-number seed even if you don’t care about reproducibility.

If you do care about reproducibility, I will mention (1) that you need to use

**sort**to put the original data in a known, reproducible order, before you generate the random variate*u*, and I will explain (2) a subtle issue that leads us to use different code for*N*≤1,000 and*N*>1,000. The code for for*N*≤1,000 is. set seed

*#*. sort*variables_that_put_data_in_unique_order*. generate double u = runiform() . sort u . keep in 1/*n*and the code for

*N*>1,000 is. set seed

*#*. sort*variables_that_put_data_in_unique_order*. generate double u1 = runiform() . generate double u2 = runiform() . sort u1 u2 . keep in 1/*n*You can use the

*N*>1,000 code for the*N*≤1,000 case. - To draw without replacement a
*P*-percent random sample, type. set seed

*#*. keep if runiform() <=*P*/100There’s no issue in this case when

*N*is large.

As I mentioned, we’ll discuss drawing random samples *with replacement* next time. Today, the topic is random samples *without replacement*. Let’s start.

**Shuffling data**

I have a deck of 52 cards, in order, the first four of which are

. list in 1/4 +-------------+ | rank suit | |-------------| 1. | Ace Club | 2. | 1 Club | 3. | 2 Club | 4. | 3 Club | +-------------+

Well, actually I just have a Stata dataset with observations corresponding to playing cards. To shuffle the deck — to place the observations in random order — type

. generate double u = runiform() . sort u

Having done that, here’s your hand,

. list in 1/5 +----------------------------+ | rank suit u | |----------------------------| 1. | Queen Club .0445188 | 2. | 5 Diamond .0580662 | 3. | 7 Club .0610638 | 4. | King Heart .0907986 | 5. | 6 Spade .0981878 | +----------------------------+

and here’s mine:

. list in 6/10 +---------------------------+ | rank suit u | |---------------------------| 6. | 8 Diamond .1024369 | 7. | 5 Club .1086679 | 8. | 8 Spade .1091783 | 9. | 2 Spade .1180158 | 10. | Ace Club .1369841 | +---------------------------+

All I did was generate random numbers — one per observation (card) — and then place the observations in ascending order of the random values. Doing that is equivalent to shuffling the deck. I used **runiform()** random numbers, meaning rectangularly distributed random numbers over [0, 1), but since I’m only exploiting the random-numbers’ ordinal properties, I could have used random numbers from any continuous distribution.

This simple, elegant, and obvious solution to shuffling data will play an important part of the solution to drawing observations without replacement. I have already more than hinted at the solution when I showed you your hand and mine.

**Drawing n observations without replacement**

Drawing without replacement is exactly the same problem as dealing cards. The solution to the physical card problem is to shuffle the cards and then draw the top cards. The solution to randomly selecting *n* from *N* observations is to put the *N* observations in random order and keep the first *n* of them.

. use cards, clear . generate double u = runiform() . sort u . keep in 1/5 (47 observations deleted) . list +---------------------------+ | rank suit u | |---------------------------| 1. | Ace Diamond .0064866 | 2. | 6 Heart .0087578 | 3. | King Spade .014819 | 4. | 3 Spade .0955155 | 5. | King Diamond .1007262 | +---------------------------+ . drop u

**Reproducibility**

You might later want to reproduce the analysis, meaning you do not want to draw another random sample, but you want to draw the *same* random sample. Perhaps you informally distributed some preliminary results and, of course, then discovered a mistake. You want to redistribute updated results and show that your mistake didn’t change results by much, and to drive home the point, you want to use the same samples as you used previously.

Part of the solution is to set the random-number seed. You might type

. set seed 49983882 . use cards, clear . generate double u = runiform() . sort u . keep in 1/5

See **help set seed** in Stata. As a quick review, when you set the random-number seed, you set Stata’s random-number generator into a fixed, reproducible state, which is to say, the sequence of random numbers that **runiform()** produces is a function of the seed. Set the seed today to the same value as yesterday, and **runiform()** will produce the same sequence of random numbers today as it did yesterday. Thus, after setting the seed, if you repeat today exactly what you did yesterday, you will obtain the same results.

So imagine that you set the random number seed today to the value you set it to yesterday and you repeat the above commands. Even so, you might not get the same results! You will not get the same results if the observations in cards.dta are in a different order yesterday and today. Setting the seed merely ensures that if yesterday the smallest value of **u** was in observation 23, it will again be in observation 23 today (and it will be the same value). If yesterday, however, observation 23 was the 6 of Clubs, and today it’s the 7 of Hearts, then today you will select the 7 of Hearts in place of the 6 of Clubs.

So make sure the data are in the same order. One way to do that is put the dataset in a known order before generating the random values on which you will sort. For instance,

. set seed 49983882 . use cards, clear . sort rank suit . generate double u = runiform() . sort u . keep in 1/5

An even better solution would add the line

. by rank suit: assert _N==1

just before the **generate**. That line would check whether sorting on variables rank and suit uniquely orders the observations.

With cards.dta, you can argue that the **assert** is unnecessary, but not because you know each rank-suit combination occurs once. You have only my assurances about that. I recommend you never trust anyone’s assurances about data. In this case, however, you can argue that the **assert** is unnecessary because we sorted on *all the variables in the dataset* and thus uniqueness is not required. Pretend there are two Ace of Clubs in the deck. Would it matter that the first card was Ace of Clubs followed by Ace of Clubs as opposed to being the other way around? Of course it would not; the two states are indistinguishable.

So let’s assume there is another variable in the dataset, say whether there was a grease spot on the back of the card. Yesterday, after sorting, the ordering might have been,

+---------------------------------+ | rank suit grease u | |---------------------------------| 1. | Ace Club yes .6012949 | 2. | Ace Club no .1859054 | +---------------------------------+

and today,

+---------------------------------+ | rank suit grease u | |---------------------------------| 1. | Ace Club no .6012949 | 2. | Ace Club yes .1859054 | +---------------------------------+

If yesterday you selected the Ace of Clubs without grease, today you would select the Ace of Clubs with grease.

My recommendation is (1) sort on whatever variables put the data into a unique order, and then verify that, or (2) sort on all the variables in the dataset and then don’t worry whether the order is unique.

**Ensuring a random ordering**

Included in our reproducible solution but omitted from our base solution was setting the random-number seed,

. set seed 49983882

Setting the seed is important even if you don’t care about reproducibility. Each time you launch Stata, Stata sets the same random-number seed, namely 123456789, and that means that **runiform()** generates the same sequence of random numbers, and that means that if you generated all your random samples right after launching Stata, you would always select the same observations, at least holding *N* constant.

So set the seed, but don’t set it too often. You set the seed once per problem. If I wanted to draw 10,000 random samples from the same data, I could code:

usedataset, clear set seed 1702213 sortvariables_that_put_data_in_unique_orderpreserve forvalues i=1(1)10000 { generate double u = runiform() sort u keep in 1/ndrop u save sample`i', replace restore, preserve }

In the example I save each sample in a file. In real life, I seldom (never) save the samples; I perform whatever analysis on the samples I need and save the results, which I usually append into a single dataset. I don’t need to save the individual samples because I can recreate them.

**And the result still might not be reproducible …**

**runiform()** draws random-numbers *with* replacement. It is thus possible that two or more observations could have the same random values associated with them. Well yes, you’re thinking, I see that it’s possible, but surely it’s so unlikely that it just doesn’t happen. But it does happen:

. clear all . set obs 100000 obs was 0, now 100000 . generate double u = runiform() . by u, sort: assert _N==1 1 contradiction in 99999 by-groups assertion is false r(9);

In the 100,000-observation dataset I just created, I got a duplicate! By the way, I didn’t have to look hard for such an example, I got it the first time I tried.

I have three things I want to tell you:

- Duplicates happen more often than you might guess.

- Do not panic about the duplicates. Because of how Stata is written, duplicates do not lower the quality of the sample selected. I’ll explain.

- Duplicates do interfere with reproducibility, however, and there is an easy way around that problem.

Let’s start with the chances of observing duplicates. I mentioned in passing last time that **runiform()** is a 32-bit random-number generator. That means **runiform()** can return any of 2^{32} values. Their values are, in order,

0 = 0 1/2^{32}= 2.32830643654e-10 2/2^{32}= 4.65661287308e-10 3/2^{32}= 6.98491930962e-10 . . . (2^{32}-2)/2^{32}= 0.9999999995343387 (2^{32}-1)/2^{32}= 0.9999999997671694

So what are the chances that in *N* draws with replacement from an urn containing these 2^{32} values, that all values are distinct? The probability *p* that all values are distinct is

2^{32}* (2^{32}-1) * ... *(2^{32}-N) p = ---------------------------- N*2^{32}

Here are some values for various values of *N*. *p* is the probability that all values are unique, and 1-*p* is the probability of observing one or more repeated values.

------------------------------------ N p 1-p ------------------------------------ 50 0.999999715 0.000000285 500 0.999970955 0.000029045 1,000 0.999883707 0.000116293 5,000 0.997094436 0.002905564 10,000 0.988427154 0.011572846 50,000 0.747490440 0.252509560 100,000 0.312187988 0.687812012 200,000 0.009498117 0.990501883 300,000 0.000028161 0.999971839 400,000 0.000000008 0.999999992 500,000 0.000000000 1.000000000 ------------------------------------

In shuffling cards we generated *N*=52 random values. The probability of a repeated values is infinitesimal. In datasets of *N*=10,000, I expect to see repeated values 1% of the time. In datasets of *N*=50,000, I expect to see repeated values 25% of the time. By *N*=100,000, I expect to see repeated values more often than not. By *N*=500,000, I expect to see repeated value in virtually all sequences.

Even so, I promised you that this problem does not affect the randomness of the ordering. It does not because of how Stata’s **sort** command is written. Remember the basic solution,

. usedataset, clear . generate double u = runiform() . sort u . keep in 1/n

Did you know **sort** has its own, private random-number generator built into it? It does, and **sort** uses its random-number generator to determine the order of tied observations. In the manuals we at StataCorp are fond of writing, “the ties will be ordered randomly” and a few sophisticated users probably took that to mean, “the ties will be ordered in a way that we at StataCorp do not know and even though they might be ordered in a way that will cause a bias in the subsequent analysis, because we don’t know, we’ll ignore the possibility.” But we meant it when wrote that the ties will be ordered randomly; we know that because we put a random number generator into **sort** to ensure the result. And that is why I can now write that repeated values of the **runiform()** function cause a reproducibility issue, but not a statistical issue.

The solution to the reproducibility issue is to draw two random numbers and use the random-number pair to order the observations:

. usedataset, clear . sortvarnames. set seed#. generate double u1 = runiform() . generate double u2 = runiform() . sort u1 u2 . keep in 1/n

You might wonder if we would ever need three random numbers. It is very unlikely. *p*, the probability of no problem, equals 1 to at least 5 digits for *N*=500,000. Of course, the chances of duplication are always nonzero. If you are concerned about this problem, you could add an **assert** to the code to verify that the two random numbers together do uniquely identify the observations:

. usedataset, clear . sortvarnames. set seed#. generate double u1 = runiform() . generate double u2 = runiform() . sort u1 u2 . by u1 u2: assert _N==1// added line. keep in 1/n

I do not believe that doing that is necessary.

**Is using doubles necessary?**

In the generation of random numbers in all of the above, note that I am storing them as doubles. For the reproducibility issue, that is important. As I mentioned in part 1, the 32-bit random numbers that **runiform()** produces will be rounded if forced into 23-bit **float**s.

Above I gave you a table of probabilities *p* that, in creating

. generate double u = runiform()

the values of **u** would be distinct. Here is what would happen if you instead stored **u** as a **float**:

u stored as --------------------------------------------------------- -------- double ---------- ----------float ---------- N p 1-p p 1-p ------------------------------------------------------------------- 50 0.999999715 0.000000285 0.999853979 0.000146021 500 0.999970955 0.000029045 0.985238383 0.014761617 1,000 0.999883707 0.000116293 0.942190868 0.057809132 5,000 0.997094436 0.002905564 0.225346930 0.774653070 10,000 0.988427154 0.011572846 0.002574145 0.997425855 50,000 0.747490440 0.252509560 0.000000000 1.000000000 100,000 0.312187988 0.687812012 0.000000000 1.000000000 200,000 0.009498117 0.990501883 0.000000000 1.000000000 300,000 0.000028161 0.999971839 0.000000000 1.000000000 400,000 0.000000008 0.999999992 0.000000000 1.000000000 500,000 0.000000000 1.000000000 0.000000000 1.000000000 -------------------------------------------------------------------

**Drawing without replacement P-percent random samples**

We have discussed drawing without replacement *n* observations from *N* observations. The number of observations selected has been fixed. Say instead we wanted to draw a 10% random sample, meaning that we independently allow each observation to have a 10% chance of appearing in our sample. In that case, the final number of observations is *expected* to be 0.1**N*, but it may (and probably will) vary from that. The basic solution for drawing a 10% random sample is

. keep if runiform() <= 0.10

and the basic solution for drawing a *P*% random sample is

. keep if runiform() <=P/100

It is unlikely to matter whether you code <= or < in the comparison. As you now know, **runiform()** produces values drawn from 2^{32} possible values, and thus the chance of equality is 2^{-32} or roughly 0.000000000232830644. If you want a *P*% sample, however, theory says you should code <=.

If you care about reproducibility, you should expand the basic solution to read,

. set seed#. usedata, clear . sortvariables_that_put_data_in_unique_order. keep if runiform() <=P/100

Below I draw a 10% sample from the card.dta:

. set seed 838 . use cards, clear . sort rank suit . keep if runiform() <= 10/100 (46 observations deleted) . list +-----------------+ | rank suit | |-----------------| 1. | 2 Diamond | 2. | 2 Heart | 3. | 3 Club | 4. | 5 Heart | 5. | Jack Diamond | |-----------------| 6. | Queen Spade | +-----------------+

**We’re not done, but we’re done for today**

In part 3 of this series I will discuss drawing random samples *with* replacement.

Pingback: The Stata Blog » Using Stata’s random-number generators, part 3, drawing with replacement()