3

我正在学习一些对书中缺失数据进行统计的技术（ Little 和 Rubin 的Statistical Analysis with Missing Data）。处理单调无响应数据的一项特别有用的功能是Sweep Operator（详细信息见第 148-151 页）。我知道 R 模块gmm具有执行此操作的swp函数，但我想知道是否有人在 Python 中实现了此函数，非常适合 Numpy 矩阵保存输入数据。我搜索了 StackOverflow，也进行了几次网络搜索，但均未成功。谢谢你的帮助。

更新：我在想我只是想念它，而不是有人会为我写它。但这也很棒。这是定义。

如果一个 PxP 对称矩阵 G 被另一个对称 PxP 矩阵 H 替换，其元素定义如下： 那么是这样，但现在遵循一个代码块：

    h_kk = -1/g_kk
    h_jk = h_kj = g_jk/g_kk for j != k
    h_jl = g_jl - g_jk g_kl / g_kk j != k, l != k
        
        
    G = [g11, g12, g13
         g12, g22, g23
         g13, g23, g33]   
    H = SWP(1,G) = [-1/g11, g12/g11, g13/g11
                   g12/g11, g22-g12^2/g11, g23-g13*g12/g11
                   g13/g11, g23-g13*g12/g11, g33-g13^2/g11]
    kvec = [k1,k2,k3]
    SWP[kvec,G] = SWP(k1,SWP(k2,SWP(k3,G)))

    Inverse function
    H = RSW(k,G)
    h_kk = -1/g_kk
    h_jk = h_kj = -g_jk/g_kk for j != k
    h_jl = g_jk g_kl / g_kk j != k, l != k    

    G == SWP(k,RSW(k,G)) == RSW(k,SWP(k,G))
4

1 回答 1

6 
def sweep(g, k):
    g = np.asarray(g)
    n = g.shape[0]
    if g.shape != (n, n):
        raise ValueError('Not a square array')
    if not np.allclose(g - g.T, 0):
        raise ValueError('Not a symmetrical array')
    if k >= n:
        raise ValueError('Not a valid row number')
    #  Fill with the general formula
    h = g - np.outer(g[:, k], g[k, :]) / g[k, k]
    # h = g - g[:, k:k+1] * g[k, :] / g[k, k]
    # Modify the k-th row and column
    h[:, k] = g[:, k] / g[k, k]
    h[k, :] = h[:, k]
    # Modify the pivot
    h[k, k] = -1 / g[k, k]
    return h

我没有办法测试上面的代码，但是我在这里找到了一个替代描述，它对非对称矩阵有效，可以计算如下：

def sweep_non_sym(a, k):
    a = np.asarray(a)
    n = a.shape[0]
    if a.shape != (n, n):
        raise ValueError('Not a square array')
    if k >= n:
        raise ValueError('Not a valid row number')
    #  Fill with the general formula
    b = a - np.outer(a[:, k], a[k, :]) / a[k, k]
    # b = a - a[:, k:k+1] * a[k, :] / a[k, k]
    # Modify the k-th row and column
    b[k, :] = a[k, :] / a[k, k]
    b[:, k] = -a[:, k] / a[k, k]
    # Modify the pivot
    b[k, k] = 1 / a[k, k]
    return b

这确实为该链接中的示例提供了正确的结果：

>>> a = [[2,4],[3,1]]
>>> sweep_non_sym(a, 0)
array([[ 0.5,  2. ],
       [-1.5, -5. ]])
>>> sweep_non_sym(sweep_non_sym(a, 0), 1)
array([[-0.1,  0.4],
       [ 0.3, -0.2]])
>>> np.dot(a, sweep_non_sym(sweep_non_sym(a, 0), 1))
array([[  1.00000000e+00,   0.00000000e+00],
       [  5.55111512e-17,   1.00000000e+00]])
于 2013-04-02T17:12:42.930 回答