1

这可能会导致 scipy/numpy，但现在我对任何功能都很满意，因为我在这些包中找不到任何东西。我有一个矩阵，其中包含多变量分布的数据（为了好玩，假设为 2）。是否有任何函数可以计算（更高）时刻？我只能找到 numpy.mean() 和 numpy.cov() ：o

谢谢 ：）

/编辑：

所以更多细节：我有多元数据，即行显示变量和列观察的矩阵。现在我想有一种简单的方法来计算该数据的联合矩，如http://en.wikipedia.org/wiki/Central_moment#Multivariate_moments中所定义。我对 python/scipy 很陌生，所以我不确定我是否是编写这个代码的最佳人选，尤其是对于 n 变量的情况（请注意，维基百科的定义是针对 n=2），而我有点期望有一些开箱即用的东西可以使用，因为我认为这将是一个标准问题。

/编辑2：

只是为了将来，如果有人想做类似的事情，下面的代码（仍在审查中）应该给出与原始矩 E(X^2)、E(Y^2) 等等效的样本。它目前仅适用于两个变量，但如果需要，它应该是可扩展的。如果您看到一些错误或不干净/unpython-nish 代码，请随时发表评论。

from numpy import *

# this function should return something as
# moments[0] = 1
# moments[1] = mean(X), mean(Y)
# moments[2] = 1/n*X'X, 1/n*X'Y, 1/n*Y'Y
# moments[3] = mean(X'X'X), mean(X'X'Y), mean(X'Y'Y),
#   mean(Y'Y'Y)
# etc
def getRawMoments(data, moment, axis=0):
    a = moment
    if (axis==0):
        n = float(data.shape[1])
        X = matrix(data[0,:]).reshape((n,1))
        Y = matrix(data[1,:]).reshape((n,1))
    else:
        n = float(data.shape[0])
        X = matrix(data[:,0]).reshape((n,1))
        Y = matrix(data[:,1]).reshape((n,11))

    result = 1
    Z = hstack((X,Y))
    iota = ones((1,n))
    moments = {}
    moments[0] = 1
    #first, generate huge-ass matrix containing all x-y combinations
    # for every power-combination k,l such that k+l = i
    # for all 0 <= i <= a
    for i in arange(1,a):
        if i==2:
            moments[i] = moments[i-1]*Z
        # if even, postmultiply with X. 
        elif i%2 == 1:
            moments[i] = kron(moments[i-1], Z.T)
        # Else, postmultiply with X.T
        elif i%2==0:
            temp = moments[i-1]
            temp2 = temp[:,0:n]*Z
            temp3 = temp[:,n:2*n]*Z
            moments[i] = hstack((temp2, temp3)) 
    # since now we have many multiple moments
    # such as x**2*y and x*y*x, filter non-distinct elements
    momentsDistinct = {}
    momentsDistinct[0] = 1
    for i in arange(1,a):
        if i%2 == 0:
            data = 1/n*moments[i]
        elif i == 1:
            temp = moments[i]
            temp2 = temp[:,0:n]*iota.T
            data = 1/n*hstack((temp2))            
        else:
            temp = moments[i]
            temp2 = temp[:,0:n]*iota.T
            temp3 = temp[:,n:2*n]*iota.T
            data = 1/n*hstack((temp2, temp3))
        momentsDistinct[i] = unique(data.flat)
    return momentsDistinct(result, axis=1)
4

0 回答 0