python中的多变量累积量和矩

Question

在Mathematica中，我可以使用MomentConvert在累积量中转换多变量矩：

MomentConvert[Cumulant[{2, 2,1}], "Moment"] // TraditionalForm

正如可以在wolframcloud中尝试的那样。

我想在python中做同样的事情。在python中有没有能够做到这一点的库？

Answer 1

至少有一个方向是我现在自己编写的：

# from http://code.activestate.com/recipes/577211-generate-the-partitions-of-a-set-by-index/

from collections import defaultdict

class Partition:

    def __init__(self, S):        
        self.data = list(S)        
        self.m = len(S)
        self.table = self.rgf_table()

    def __getitem__(self, i):
        #generates set partitions by index
        if i > len(self) - 1:
             raise IndexError
        L =  self.unrank_rgf(i)
        result = self.as_set_partition(L)
        return result

    def __len__(self):
        return self.table[self.m,0]

    def as_set_partition(self, L):
        # Transform a restricted growth function into a partition
        n = max(L[1:]+[1])
        m = self.m
        data = self.data
        P = [[] for _ in range(n)]
        for i in range(m):
            P[L[i+1]-1].append(data[i])
        return P

    def rgf_table(self):
        # Compute the table values 
        m = self.m
        D = defaultdict(lambda:1)
        for i in range(1,m+1):
            for j in range(0,m-i+1):
                D[i,j] = j * D[i-1,j] + D[i-1,j+1]
        return D

    def unrank_rgf(self, r):
        # Unrank a restricted growth function
        m = self.m
        L = [1 for _ in range(m+1)]
        j = 1
        D = self.table
        for i in range(2,m+1):
            v = D[m-i,j]
            cr = j*v
            if cr <= r:
                L[i] = j + 1
                r -= cr
                j += 1
            else:
                L[i] = r // v + 1
                r  %= v
        return L



# S = set(range(4))
# P = Partition(S)
# for x in P:
#      print (x)


# using https://en.wikipedia.org/wiki/Cumulant#Joint_cumulants 
import math

def Cum2Mom(arr, state):
    def E(op):
        return qu.expect(op, state) 

    def Arr2str(arr,sep):
        r = ''
        for i,x in enumerate(arr):
            r += str(x)
            if i<len(arr)-1:
                r += sep                
        return r

    if isinstance( arr[0],str):
        myprod = lambda x: Arr2str(x,'*')
        mysum = lambda x: Arr2str(x,'+')
        E=lambda x: 'E('+str(x)+')'
        myfloat = str
    else:
        myfloat = lambda x: x
        myprod = np.prod
        mysum = sum
    S = set(range(len(arr)))
    P = Partition(S)  

    return mysum([  
        myprod([myfloat(math.factorial(len(pi)-1)   * (-1)**(len(pi)-1))
        ,myprod([
            E(myprod([
                arr[i]
                for i in B
            ]))
            for B in pi])])
        for pi in P])

print(Cum2Mom(['a','b','c','d'],1)    )
import qutip as qu
print(Cum2Mom([qu.qeye(3) for i in range(3)],qu.qeye(3))    )

它被设计为使用qutip对象，也可以使用字符串来验证正确的分隔和前置因子。

变量的指数可以通过重复变量来表示。