3

本から欠損データを使用して統計を行うためのいくつかの手法を学んでいます（リトルとルービンによる欠損データを使用した統計分析）。単調な無応答データを操作するのに特に役立つ機能の 1 つは、Sweep Operator (詳細は 148 ～ 151 ページ) です。Rモジュールgmmにはこれを行うswp関数があることは知っていますが、理想的にはNumpy行列が入力データを保持するために、Pythonでこの関数を実装した人がいるかどうか疑問に思っていました。私は StackOverflow を検索し、いくつかの Web 検索も成功しませんでした。助けてくれてありがとう。

更新: 誰かが私のために書いてくれるとは思っていませんでした。しかし、それも素晴らしいでしょう。これが定義です。

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

上記のコードをテストする方法はありませんが、次のように計算できる非対称行列に有効な別の説明hereを見つけました。

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 に答える