基于递归C（n，k）= C（n - 1，k）+ C（n - 1，k - 1），记忆化，使用numpy操作。

import numpy as np

def binnings(n, k, cache={}):
    if n == 0:
        return np.zeros((1, k))
    if k == 0:
        return np.empty((0, 0))
    args = (n, k)
    if args in cache:
        return cache[args]
    a = binnings(n - 1, k, cache)
    a1 = a + (np.arange(k) == 0)
    b = binnings(n, k - 1, cache)
    b1 = np.hstack((np.zeros((b.shape[0], 1)), b))
    result = np.vstack((a1, b1))
    cache[args] = result
    return result

if __name__ == '__main__':
    import timeit
    print timeit.timeit('binnings(20, 5, {})', setup='from __main__ import binnings', number=1)

字符串
（20，5）的时间（秒）：

0.0129251480103

型

在numpy中使用一些不同的技巧可能会有更快的方法。numpy.indices是你想要开始的地方。一旦你将它与rollaxis合并结合起来，它本质上相当于itertools.product。Sven Marnach的答案in this question是一个很好的例子（在他的最后一个例子中有一个小错误，但是，这是你想要使用的。它应该是numpy.indices((len(some_list) + 1), * some_length...）
然而，对于这样的东西，使用itertools可能会更具可读性。
numpy.fromiter比显式转换为列表要快一些，特别是当你给予迭代器中元素的数量时。主要的优点是它使用的内存少得多，这在大迭代器的情况下非常有用。
下面是一些使用numpy.indices技巧和各种将迭代器转换为numpy数组的方法的例子：

import itertools
import numpy as np
import scipy.special

def fixed_total_product(bins, num_items):
    return itertools.ifilter(lambda combo: sum(combo) == num_items,
            itertools.product(xrange(num_items + 1), repeat=bins))

def fixed_total_product_fromiter(bins, num_items):
    size = scipy.special.binom(bins - 1 + num_items, num_items)
    combinations = fixed_total_product(bins, num_items)
    indicies = (idx for row in combinations for idx in row)
    arr = np.fromiter(indicies, count=bins * int(size), dtype=np.int)
    return arr.reshape(-1, bins)

def fixed_total_product_fromlist(bins, num_items):
    return np.array(list(fixed_total_product(bins, num_items)), dtype=np.int)

def fixed_total_product_numpy(bins, num_items):
    arr = np.rollaxis(np.indices((num_items + 1,) * bins), 0, bins + 1)
    arr = arr.reshape(-1, bins)
    arr = np.arange(num_items + 1)[arr]
    sums = arr.sum(axis=1)
    return arr[sums == num_items]

m, n = 5, 20

if __name__ == '__main__':
    import timeit
    list_time = timeit.timeit('fixed_total_product_fromlist(m, n)',
            setup='from __main__ import fixed_total_product_fromlist, m, n',
            number=1)
    iter_time = timeit.timeit('fixed_total_product_fromiter(m, n)',
            setup='from __main__ import fixed_total_product_fromiter, m, n',
            number=1)
    numpy_time = timeit.timeit('fixed_total_product_numpy(m, n)',
            setup='from __main__ import fixed_total_product_numpy, m, n',
            number=1)

    print 'All combinations of {0} items drawn from a set of {1} items...'.format(m,n)
    print '  Converting to a list and then an array: {0} sec'.format(list_time)
    print '  Using fromiter: {0} sec'.format(iter_time)
    print '  Using numpy.indices: {0} sec'.format(numpy_time)

字符串
至于时间：

All combinations of 5 items drawn from a set of 20 items...
  Converting to a list and then an array: 2.75901389122 sec
  Using fromiter: 2.10619592667 sec
  Using numpy.indices: 1.44955015182 sec

型
你会注意到它们都不是特别快。
您使用的是蛮力算法（生成所有可能的组合，然后过滤它们），我只是在这里的基于numpy的示例中复制它。
可能有一个更有效的方法来做到这一点！然而，我不是一个CS或数学的人无论如何，所以我不知道是否有一个众所周知的算法来做到这一点，而不是首先生成所有可能的组合.
祝你好运，无论如何！

我知道我在这里恢复了一个很老的线程，但我希望一个好的解决方案仍然会受到赞赏（即使不再是OP）。
这个问题本身类似于在代码中表示多变量多项式。例如（对应于3个bin中的4个项目的例子），你从扩展(x+y+z)^4得到的多项式类似于x^4*y^0*z^0 + 4*x^3*y^1*z^0 + ...，x，y和z的指数是数组的元素：

[4,0,0],
[3,1,0],
[3,0,1],
[2,2,0],
[2,1,1],
[2,0,2],
[1,3,0],
[

字符串
如果你仔细观察，你会看到例如3右边的元素是(y+z)^1的指数，2右边的元素是(y+z)^2的指数等等。更一般地，在p右边你有(y+z)^(4-p)的指数。这表明了一些递归结构：假设e（m，p）表示m变量的指数（= bin）和p总指数（=items），然后通过从0到p中选择第一个指数q来获得所有指数，剩余的m-1变量的指数由e(m-1,p-q)给出。
在python和numpy中，您可以像这样公式化：

def multiindex_exact(m, p):
    if m == 0:
        return np.zeros((1 if p==0 else 0, 0), np.int8)
    else:
        I = np.zeros((0, m), np.int8)
        for q in reversed(range(0, p + 1)):
            J = multiindex_exact(m - 1, p - q)
            Jn = np.full((J.shape[0], 1), q)
            I = np.vstack((I, np.hstack((Jn, J))))
        return I

型
当然，您可以通过预先计算数组大小并直接填充值来提高效率。
如果您需要将p项分配到m bin中，而不是精确地分配，那么您可以使用以下类似的代码。

def multiindex_lower(m, p):
    if m == 0:
        return np.zeros((1, 0), np.int8)
    else:
        I = np.zeros((0, m), np.int8)
        for q in range(0, p + 1):
            J = multiindex_lower(m - 1, q)
            Jn = q - J.sum(1).reshape((J.shape[0], 1))
            I = np.vstack((I, np.hstack((J, Jn))))
        return I

型