统计Sweep算子的Python实现

我正在学习一些用书中缺失的数据进行统计的技术（缺失数据的统计分析作者：利特尔和鲁宾）。对于处理单调无响应数据来说，一个特别有用的函数是扫频操作员（详情见第 148-151 页）。我知道 R 模块gmm有swp函数可以做到这一点，但我想知道是否有人在 Python 中实现了这个函数，理想情况下 Numpy 矩阵可以保存输入数据。我搜索了 StackOverflow，也进行了几次网络搜索，但都没有成功。谢谢你的帮助。

这是定义。

如果将 PxP 对称矩阵 G 替换为另一个对称 PxP 矩阵 H（其元素定义如下），则称其在行和列 k 上进行扫描：

    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))

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

我无法测试上面的代码，但我找到了替代描述here http://www-rci.rutgers.edu/~dhjones/APPLIED_LINEAR_STATISTICAL_MODELS%28PHD%29/LECTURES/LECTURE06/3-The%20sweep%20operator.pdf，对于非对称矩阵有效，计算公式如下：

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]])