10
$\begingroup$

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.

$\endgroup$
0
10
$\begingroup$

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)

Frequency distribution

$\endgroup$
1
  • 1
    $\begingroup$ 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). $\endgroup$ – whuber Sep 19 '12 at 14:14
7
$\begingroup$

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.

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

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.