私信  •  关注

Paul Panzer

Paul Panzer 最近创建的主题
Paul Panzer 最近回复了

这是一个基于年轻图表的方法。例如,4个篮子,6个鸡蛋,每个篮子最多3个鸡蛋。如果我们根据篮子的装满程度来订购篮子,我们就会得到年轻的图表。

x x x x   x x x x   x x x     x x x      x x
x x       x         x x x     x x        x x
          x                   x          x x

下面的代码列举了所有可能的年轻图表,并为每个图表列举了所有可能的排列。

同样的逻辑也可以用来计数。

from itertools import product, combinations
from functools import lru_cache
import numpy as np

def enum_ord_part(h, w, n, o=0):
    if h == 1:
        d = n
        for idx in combinations(range(w), d):
            idx = np.array(idx, int)
            out = np.full(w, o)
            out[idx] = o+1
            yield out
    else:
        for d in range((n-1)//h+1, min(w, n) + 1):
            for idx, higher in product(combinations(range(w), d),
                                       enum_ord_part(h-1, d, n-d, o+1)):
                idx = np.array(idx)
                out = np.full(w, o)
                out[idx] = higher
                yield out

def bc(n, k):
    if 2*k > n:
        k = n-k
    return np.prod(np.arange(n-k+1, n+1, dtype='O')) // np.prod(np.arange(1, k+1, dtype='O'))

@lru_cache(None)
def count_ord_part(h, w, n):
    if h == 1:
        return bc(w, n)
    else:
        return sum(bc(w, d) * count_ord_part(h-1, d, n-d)
                   for d in range((n-1)//h+1, min(w, n) + 1))

几个例子:

>>> for i, l in enumerate(enum_ord_part(3, 4, 6), 1):
...     print(l, end=' ' if i % 8 else '\n')
... 
[3 3 0 0] [3 0 3 0] [3 0 0 3] [0 3 3 0] [0 3 0 3] [0 0 3 3] [3 2 1 0] [2 3 1 0]
[3 1 2 0] [2 1 3 0] [1 3 2 0] [1 2 3 0] [2 2 2 0] [3 2 0 1] [2 3 0 1] [3 1 0 2]
[2 1 0 3] [1 3 0 2] [1 2 0 3] [2 2 0 2] [3 0 2 1] [2 0 3 1] [3 0 1 2] [2 0 1 3]
[1 0 3 2] [1 0 2 3] [2 0 2 2] [0 3 2 1] [0 2 3 1] [0 3 1 2] [0 2 1 3] [0 1 3 2]
[0 1 2 3] [0 2 2 2] [3 1 1 1] [1 3 1 1] [1 1 3 1] [1 1 1 3] [2 2 1 1] [2 1 2 1]
[2 1 1 2] [1 2 2 1] [1 2 1 2] [1 1 2 2]
>>> 
>>> print(f'{count_ord_part(2, 26, 30):,}')
154,135,675,070
>>> print(f'{count_ord_part(50, 30, 1000):,}')
63,731,848,167,716,295,344,627,252,024,129,873,636,437,590,711
5 年前
回复了 Paul Panzer 创建的主题 » 复数numpy数组的python均值漂移聚类

在注释/聊天中,我们至少发现了一个问题,即

(cov_w + I)^-1 @ cov_b                       (1)

不是它应该的实数,而是返回重要的虚分量。这里@是矩阵乘法,cov_w和cov_b是协方差矩阵,i是恒等矩阵。这可以通过计算(cov_w+i)^-1的矩阵平方根来解决,我们称它为sq,然后使用(1)类似于

SQ @ cov_b @ SQ                              (2)

因此具有相同的特征值,如果v是(2)的特征向量,那么(1)的(右)特征向量是sq@v。

我们得到的结果是,因为(2)是一个对称矩阵,它的特征分解可以用 numpy.linalg.eigh 这保证了纯粹真实的结果。 eigh 也可用于计算sq,请参见 here . 一定要绕过逆时针方向 哎哟 直接打开 cov_w + I 甚至 cov_w .