# How do I generate numbers according to a Soliton distribution?

The Soliton distribution is a discrete probability distribution over a set $\{1,\dots, N\}$ with the probability mass function

$$p(1)=\frac{1}{N},\qquad p(k)=\frac{1}{k(k-1)}\quad\text{for }k\in\{2,\dots, N\}$$

I'd like to use it as part of an implementation of an LT code, ideally in Python where a uniform random number generator is available.

If we start at $k=2$, the sums telescope, giving $1-1/k$ for the (modified) CDF. Inverting this, and taking care of the special case $k=1$, gives the following algorithm (coded in R, I'm afraid, but you can take it as pseudocode for a Python implementation):

rsoliton <- function(n.values, n=2) {
x <- runif(n.values)         # Uniform values in [0,1)
i <- ceiling(1/x)            # Modified soliton distribution
i[i > n] <- 1                # Convert extreme values to 1
i
}


As an example of its use (and a test), let's draw $10^5$ values for $N=10$:

n.trials <- 10^5
i <- rsoliton(n.trials, n=10)
freq <- table(i) / n.trials  # Tabulate frequencies
plot(freq, type="h", lwd=6)


• For the related "robust" soliton distribution, you probably have to settle for a slightly less efficient solution (based on a binary search or the equivalent). – whuber Sep 19 '12 at 14:14

# Python (adapted from @whuber's R solution)

from __future__ import print_function, division
import random
from math import ceil

def soliton(N, seed):
prng = random.Random()
prng.seed(seed)
while 1:
x = random.random() # Uniform values in [0, 1)
i = int(ceil(1/x))       # Modified soliton distribution
yield i if i <= N else 1 # Correct extreme values to 1

if __name__ == '__main__':
N = 10
T = 10 ** 5 # Number of trials
s = soliton(N, s = soliton(N, random.randint(0, 2 ** 32 - 1)) # soliton generator
f = [0]*N                       # frequency counter
for j in range(T):
i = next(s)
f[i-1] += 1

print("k\tFreq.\tExpected Prob\tObserved Prob\n");

print("{:d}\t{:d}\t{:f}\t{:f}".format(1, f[0], 1/N, f[0]/T))
for k in range(2, N+1):
print("{:d}\t{:d}\t{:f}\t{:f}".format(k, f[k-1], 1/(k*(k-1)), f[k-1]/T))


## Sample Output

k   Freq.   Expected Prob   Observed Prob

1   9965    0.100000    0.099650
2   49901   0.500000    0.499010
3   16709   0.166667    0.167090
4   8382    0.083333    0.083820
5   4971    0.050000    0.049710
6   3354    0.033333    0.033540
7   2462    0.023810    0.024620
8   1755    0.017857    0.017550
9   1363    0.013889    0.013630
10  1138    0.011111    0.011380


## Requirements

The code should work in Python 2 or 3.

• @whuber LT implementation now on GitHub. Not perfect, but it's a start. – Alex Chamberlain Sep 20 '12 at 14:16