I'm trying to run some basic simulations of if one thing can randomly guess another (is the intersection of two sets non-null) and I'm curious if different implementations in software (which are very different in terms of performance) are mathematically identical:

  1. Generate two random sets of numbers (from among $n$ possible values), of sizes $x$ and $y$. Find the intersection. Is it null?
  2. Generate one random set of numbers, size $x$. Randomly pick a consecutive range of values, $y$ long. Do any values in the set fall within that random range?
  3. Generate one random set of numbers, size $x$. Do any of the values fall within the first $y$ possible values?
  • $\begingroup$ I assume you mean to choose your numbers from a uniform distribution on the total interval. Given that, then your three problems are the same. There is nothing to distinguish any one collection of $y$ elements from any other. Of course, if the distribution is non-uniform then this is false. $\endgroup$ – lulu Aug 22 '15 at 3:20

There are often three issues with a simulation.

First, is the pseudoranom generator (PRNG) adequate. (a) Does it produce values that are not distinguishable from random for all practical purposes. The numbers should be independent. Also, most PRNGs give numbers that are supposed to be uniformly distributed in the interval (0, 1). For example, each interval of length 1/10 should get about 1/10 of the values.

Most modern statistical software uses high quality PRNGs. For example, R uses the Mersenne-twister, which has been vetted for good behavior through an impressive battery of benchmark tests. If you have R, you can read more about PRNGs by typing ? .Random.seed into the Console window. (Some older software still retains inadequate generator.)

Each simulation can gives a different result. That is in the nature of randomness. An exception is that you can specify a seed in a PRNG; then the same seed, the same PRNG, and the same program will always give the same result.

Second, you need to be able to program the particular random experiment you have in mind. The software should have a way to convert random values in (0, 1) to whatever kind of probability model you have in mind. Most statistical software has a lot of pre-programmed distributions. If you want to use one of the more esoteric distributions, you might have to do some of this programming for yourself.

Different software packages can use different methods to program the same model. Proper programming can lead to the same statistical model in any software, subject to the capability of software and hardware.

Third, often you want to perform a random experiment many thousands of times. Then it is essential to be able to determine results from each run automatically and to store results so you can summarize them later.

Now to discuss the first of the specific probability experiments you mentioned in your question as an example of these issues. (I have not seen this problem before.)

Suppose I roll fair dice, so that $n = 6.$ Maybe one of them $x = 2$ times and the other $y = 3$ times. Call the two sets $a$ and $b$. Then in R, you could easily do this experiment once as follows:

 a = sample(1:6, 2, repl=T);  b = sample(1:6, 4, repl=T)
 a;  b
 ##  4 5
 ##  2 5 6 3

It is easy to see that the intersection is not null. But if I am running 100,000 of these experiments to see the proportion of time the intersection is null, then I need an $automated$ way to test that. Here is one way:

 a = sample(1:6, 2, repl=T);  b = sample(1:6, 4, repl=T)
 a; b
 ##  1 4
 ##  3 1 4 3
 match(a, b, nomatch=0)
 ## 2 3                   # indexes in b matching a
 match(a, b, nomatch=0) > 0
 ##  TRUE TRUE            # TRUE indicates match
 sum(match(a, b, nomatch=0) > 0)
 ##  2                    # sum of logical vector counts TRUE's

So all that's really necessary is:

 a = sample(1:6, 2, repl=T);  b = sample(1:6, 4, repl=T)
 sum(match(a, b, nomatch=0) > 0)
 ## 0
 a = sample(1:6, 2, repl=T);  b = sample(1:6, 4, repl=T)
 sum(match(a, b, nomatch=0) > 0)
 ## 2
 a = sample(1:6, 2, repl=T);  b = sample(1:6, 4, repl=T)
 sum(match(a, b, nomatch=0) > 0)
 ## 1

Now I write a loop to store results from many runs.

 m = 100000;  a.in.b = numeric(m)
 for (i in 1:m) {     
   a = sample(1:6, 2, repl=T);  b = sample(1:6, 4, repl=T)
   a.in.b[i] = sum(match(a, b, nomatch=0) > 0) }
 mean(a.in.b == 0)   # proportion mutually exclusive
 ##  0.24425

About a quarter of the cases yield mutually exclusive sequences. With m = 10^6 runs of the experiment, I could get the probability of 'mutually exclusive' to about 3 places.

Note: I think your proposed models 2 and 3 need a little more attention before they can be turned directly into specific experiments that are clearly different from model 1.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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