Deterministic sampling from discrete distribution

Question

I'm working on a generalization of the Min-Hash algorithm to allow the meaningful comparison of ordered values such as integers. The core trick is to use deterministic randomness as a replacement for hash-functions.

However, to get this thing to work I need a way to deterministically sample from a discrete distribution. I.e. if I repeatedly draw samples from a discrete distribution I need to get the same set of values for every call.

Library functions like sample of the languageR use standard sampling algorithms which typically use uniform-sampling from the cumulative distribution. Hence, they need some pseudo-random numbers. By seeding and resetting this random number generator prior to each call, it is possible to cut out the randomness and make things deterministic. sample then behaves like a mathematical function.

Consider following example:

value:    a       b      c     d    e     f     g      h     i    j
prob:    0.070 0.0774 0.083 0.090 0.096 0.103 0.109 0.116 0.122 0.129

setSeed(1)
sample(value,5, prob, replace = TRUE)
>>> a, c, e, j, a
setSeed(1)
sample(value,5, prob, replace = TRUE)
>>> a, c, e, j, a

However, in case the distribution changes slightly (i.e another value k gets inserted) the values returned from these standard sampling algorithms change drastically.

setSeed(1)
sample(value, 5, prob, replace = TRUE)
>>> a, c, e, j, a
insert('k')
>>> value:    a     b     c     d    e     f     g      h     i     j    k
>>> prob:  0.062 0.068 0.073 0.079 0.085 0.090 0.096 0.102 0.107 0.113 0.119
setSeed(1)
sample(value,5, prob, replace = TRUE)
>>> b, d, f, k, b

So, despite the distribution changed only marginally, the drawn sample changed broadly. However, I would like to see just a minor change, something like

>>> a, c, e, j, a   <-- initial sample from initial density
>>> a, c, k, j, a   <-- desired sample after small density change
>>> b, d, f, k, b   <-- actual sample after small density change

A slow but working algorithm would be to 'simulate' this sampling by using 10.000 bins. Each bin contains a single 'value'. The number of bins of a certain value corresponds to the probability of the value in the distribution. Drawing from this simulation then works by drawing integers i 1-10.000 and returning the value of the bin with number i. If bins get replaced with another value, the drawn sample changes only slightly, because on every call the same bins get selected.

So the problems with standard algorithms is, that they typically sort/rearrange these bins to get the speedup. However, thats not possible in that case.

Is there a way to sample from the density distribution itself while ensuring that the change in the drawn sample is similar to the change of the distribution?

Answer 1

One general way to generate similar sequences for distributions with similar probabilities is the following. Suppose $A$ is an ordered finite alphabet $(a,b,c,\ldots)$ with a probability distribution $p_A$. To draw a value at random from $A$, generate a vector of independent uniform variates $\mathbf{U}=(U_a, U_b, U_c, \ldots)$. If $U_a \le p_A(a)$, choose $a$. Otherwise, recursively draw a value from the remaining letters $A^\prime = A-\{a\} = (b,c,\ldots)$ using the vector $\mathbf{U}^\prime=(U_b, U_c, \ldots)$ and the probabilities $$p_{A^\prime}(b) = \frac{p_A(b)}{1-p_A(a)},\ p_{A^\prime}(c) = \frac{p_A(c)}{1-p_A(a)},$$ etc.

You can re-use the same vector $\mathbf{U}$ for any other distribution on $A$.

Using this method, the expected frequency with which the same letter would be drawn from distributions $p_A$ and $q_A$ is the frequency with which $a$ would be drawn from both distributions, equal to the smaller of $p_A(a)$ and $q_A(a)$, plus the expected frequency with which the same letter would be drawn from $(b,c,\ldots)$, conditional on $a$ not being drawn from either distribution.

This method is the best you can do by assigning each letter to its own connected interval. With additional work it's possible to make the two sequences agree even more frequently, but you would have to assign the extra letter "k" to a complicated subset of $(0,1]$.

---

Here is R code to generate n symbols from an alphabet with probability vector prob.

s <- function(n, prob) {
  k <- length(prob)
  q <- prob / rev(cumsum(rev(prob)))
  u <- matrix(runif(n*k), nrow=k, byrow=TRUE) < q
  apply(u, 2, function(x) match(TRUE, x))
}

Let's generate samples of size 10,000 from distributions like those in the question. The output shows the first 60 draws from each, using the same starting seed. They are remarkably similar.

prob <- c(75, 77, 83, 90, 96, 103, 109, 116, 122, 129, rep(0, 16))
prob.k <- prob
prob.k[11] <- 0.119/(1-0.119)*sum(prob)

seed <- 17
N <- 1e4
set.seed(seed); x <- letters[s(N, prob)]
set.seed(seed); x.k <- letters[s(N, prob.k)]

rbind(First=paste0(head(x, 60), collapse=""),
      Second=paste0(head(x.k, 60), collapse=""))

Here it is:

First  "geifgefhiifafbhfhijgfcdiebjfjegajgggidghchhjjfgdjheicbhbjica"
Second "geifgefhiikafbhgiikifcdiibjfkjkakiggkdghchijjfgdkhhkkbhbkkcj"

You may check that the actual frequencies are close to the intended ones:

rbind(First=c(table(x), k=0), Second=table(x.k))

This output is

         a   b   c   d   e    f    g    h    i    j    k
First  755 775 808 842 995 1068 1111 1184 1206 1256    0
Second 657 666 693 739 872  955  973 1056 1074 1144 1171

The degree of similarity (that is, proportion of time the two sequences are expected to agree) is readily computed recursively.

similarity <- function(x, y) {
  if (min(length(x), length(y)) == 0) return (0)
  a <- min(x[1], y[1])
  x.s <- sum(x[-1])
  y.s <- sum(y[-1])
  if (x.s > 0 & y.s > 0) {
    b <- max(x[1], y[1])
    x <- x[-1]/sum(x[-1])
    y <- y[-1]/sum(y[-1])
    b <- (1-b) * similarity(x, y)
  } else {
    b <- 0
  }
  return (a + b)
}
similarity(prob/sum(prob), prob.k/sum(prob.k))

The output is

0.7568941

In fact, in this simulation the observed frequency was close to that:

mean(x.k == x)

[1] 0.754