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))
------------------------------------------------------
--
http://mail.python.org/mailman/listinfo/python-list
