My Python rewrite of Brendan Frey's product form sampling code.

pfball.py

import numpy as np

def sample_pf_ball(p,b):
    """
    Samples N binary variables using the product form distribution given in
    the vector b (of length N), where b[n] is the probability that s[n]=1.
    However, only K variables may on at a time, where K is the length of
    vector b minus 1, and b[k] is a distribution that biases the number of
    variables set to 1.
    %
    Formally, this function samples the binary vector s from
    %
        P(s)= (prod_n p(n)^s(n)(1-p(n))^(1-s(n)))*b(1+sum_n s(n))/Z,
    %
    where Z is the normalizing factor.
    %
    Copyright 2008 Brendan J Frey. Translated to Python by David Warde-Farley.
    Use freely for noncommercial purposes.
    """
    N=len(p); K=len(b)-1

    # Run backward algorithm, computing messages
    m = np.zeros((K+1,N)); m[:,N-1]=b
    for n in xrange(N-2,0,-1):
        m[K,n] = m[K,n+1]*(1-p[n+1])
        m[0:K,n] = m[0:K,n+1]*(1-p[n+1])+m[1:(K+1),n+1]
        m[:,n] /= np.sum(m[:,n])

    # Run forward algorithm, computing samples
    s = np.zeros((N,))
    j = 1 # Number of active units so far plus 1
    for n in xrange(N):
        if j-1 < K:
            s[n] = np.random.rand()<1./(1.+m[j-1,n]*(1-p[n])/m[j,n]/p[n])
            j += s[n]
        else:
            break
    return s
~