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

sample(value,5, prob, replace = TRUE)
>>> a, c, e, j, a
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.

sample(value, 5, prob, replace = TRUE)
>>> a, c, e, j, a
>>> 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
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?

  • 1
    $\begingroup$ sample samples values (pseudo-) randomly, so you should expect the output values to be rather random than predictable... sample does not use any CDF, it just returns random integers. Also: what do you mean by sampling from density? $\endgroup$
    – Tim
    Feb 16, 2016 at 13:43
  • 1
    $\begingroup$ Sorry but your code and your question are unclear. sample works in a very simple manner: it samples values (pseudo-)randomly, by default each value with the same probability, but if you use prob parameter it uses different probabilities as defined by user. Change in output is the expected behavior of this function if you change the input. Refer to the function documentation if it is not clear. $\endgroup$
    – Tim
    Feb 16, 2016 at 14:22
  • $\begingroup$ WHAT property? You have not defined what is the expected behavior of the algorithm that you are describing. $\endgroup$
    – Tim
    Feb 16, 2016 at 14:39
  • $\begingroup$ I updated the post $\endgroup$ Feb 16, 2016 at 15:55
  • 1
    $\begingroup$ I think the question remains unclear. I have the feeling that the point you are leaning to is "why don't small changes to the sampling frame result in small changes in the sample drawn, if we use the same random numbers for the selection?" But if you change the $n$ that you are drawing from, then clearly some change is inevitable. Perhaps what you really need to do is clarify what kind of change is concerning you? $\endgroup$
    – Silverfish
    Feb 16, 2016 at 16:57

1 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


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

mean(x.k == x)

[1] 0.754   
  • $\begingroup$ This approach looks really promising. I'll give it a try. Thanks a lot! $\endgroup$ Feb 16, 2016 at 19:03
  • $\begingroup$ After testing it I stumbled upon a major flaw. If new values at the left side of the distribution (i.e. at the left of the a) are inserted, the change of the sample is again broadly. Just like in other algorithms the problem is the uniform sampling from the CDF. This flaw can be omitted, by knowing the whole alphabet of values beforehand and setting their probability to 0 . However, that's not straight forward in my case. $\endgroup$ Feb 17, 2016 at 11:20
  • $\begingroup$ TLDR: removing the a results in a big change, but setting the probability to 0 is fine. $\endgroup$ Feb 17, 2016 at 11:43

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.