number generator

[Duncan Smith]

greg wrote:
> Gabriel Genellina wrote:
> 
>> The 5th number is not "random".
> 
> 
> More precisely, the fifth number is not *independent*
> of the others. You can't have five independent random
> numbers that sum to 50; only four independent numbers
> plus a dependent one.
> 
> -- 
> Greg

In the interests of precision, that should probably read "are
constrained to sum to 50"; it's quite possible to generate a sequence of
independent random variates that just happen to sum to 50 (as I'm sure
you know).

A fairly efficient way of generating random multinomial variates (which
would satisfy the OP's problem) is based on generating sequences of
independent Poisson variates (with suitably chosen parameters) and then
incrementing the counts in an appropriate fashion.  A fair chunk of code
is pasted below.  (If anyone knows of any more efficient methods I'd be
interested.)

Duncan

----------------------------------------------------

import math
import random

def multinomial_variate(probs, n):
    """
    Returns a random multinomial variate
    with parameters I{probs} and I{n}.

    @type  probs: sequence of C{floats}
    @param probs: multinomial probabilities
    @type      n: C{int}
    @param     n: total
    @rtype:       C{list} of C{floats}
    @return:      multinomial variate
    """
    k = len(probs)
    root_n = round(n**0.5)
    if k > root_n:
        lambda_hat = n - root_n - (k - root_n)**0.5
    else:
        lambda_hat = n - k
    gens = [gen_poisson(p * lambda_hat) for p in probs]
    gen_obs = gen_discrete(probs)
    while True:
        samp = [gen.next() for gen in gens]
        tot = sum(samp)
        if tot == n:
            break
        elif tot < n:
            for _ in xrange(n - tot):
                samp[gen_obs.next()] += 1
            break
    return samp

def _alias_table(probs):
    """
    Returns lists of aliases and cutoffs (see Kronmal and
    Peterson, 1979, I{The American Statistician}, Vol.33,
    No.4., pp. 214-218).

    @type  probs: sequence of I{floats}
    @param probs: sequence of probabilities
    @rtype:       C{tuple} of C{lists}
    @return:      a C{list} of aliases and a C{list}
                  of cutoffs
    """
    n = len(probs)
    r = [n * p for p in probs] # mean of r is 1
    a = [0] * n
    G = set()
    S = set()
    # partition r into values lower and higher than mean
    for i, val in enumerate(r):
        if val < 1:
            S.add(i)
        else:
            G.add(i)
    while True:
        try:
            j, k = S.pop(), G.pop()
        except KeyError:
            break
        a[j] = k
        r[k] += r[j] - 1
        if r[k] < 1:
            S.add(k)
        else:
            G.add(k)
    return a, r

def gen_discrete(probs):
    """
    Generates discrete (0, 1, ..., I{n}-1) variates from the
    density defined by I{probs} by the alias method.

    @type   probs: sequence of C{floats}
    @param  probs: probabilities of corresponding values
    @rtype:        generator of I{ints}
    @return:       generator of independent discrete variates
    """
    a, r = _alias_table(probs)
    n = len(probs)
    rand = random.random
    while True:
        rnd_n = rand() * n
        i = int(rnd_n)
        if rnd_n - i < r[i]:
            yield i
        else:
            yield a[i]

def gen_poisson(lambd):
    """
    Generates random Poisson variates with parameter I{lambd}.

    @type   lambd: C{float}
    @param  lambd: mean of Poisson pmf
    @rtype:        generator of I{ints}
    @return:       generator of independent Poisson variates
    """
    if lambd < 10:
        rand = random.random
        L = math.exp(-lambd)
        while True:
            k = 0
            p = 1
            while p >= L:
                k += 1
                p *= rand()
            yield k - 1
    else:
        alpha = int(7 * lambd / 8)
        gamma_var = random.gammavariate
        while True:
            X = gamma_var(alpha, 1)
            if X < lambd:
                new_lambd = lambd - X
                yield alpha + poisson_variate(new_lambd)
            else:
                yield binomial_variate(alpha - 1, lambd / X)

def poisson_variate(lambd):
    """
    Returns a random Poisson variate with parameter I{lambd}.

    @type   lambd: C{float}
    @param  lambd: mean of Poisson pmf
    @rtype:        I{int}
    @return:       Poisson variate
    """
    if lambd < 10:
        rand = random.random
        L = math.exp(-lambd)
        k = 0
        p = 1
        while p >= L:
            k += 1
            p *= rand()
        return k - 1
    else:
        a = int(7 * lambd / 8)
        X = random.gammavariate(a, 1)
        if X < lambd:
            new_lambd = lambd - X
            return a + poisson_variate(new_lambd)
        else:
            return binomial_variate(a - 1, lambd / X)

def binomial_variate(n, p):
    """
    Returns a random binomial variate with parameters I{n}
    and I{p}.

    @type  n: C{int}
    @param n: number of trials
    @type  p: C{float}
    @param p: probability of success
    @rtype:   I{int}
    @return:  number of successful trials
    """
    if n < 10:
        rand = random.random
        res = 0
        for _ in range(n):
            if rand() < p:
                res += 1
        return res
    else:
        a = 1 + n // 2
        b = n - a + 1
        X = random.betavariate(a, b)
        if X >= p:
            return binomial_variate(a - 1, p / X)
        else:
            return a + binomial_variate(b - 1, (p - X) / (1 - X))

------------------------------------------------------