2017-03-13 18 views
5

サイズ= nのサンプルがあります。NumPy:累積中央値を計算する

私はそれぞれiについて計算したい:1 < = i < = n numpyのsample[:i]の中央値。 は例えば、私は、私はそれぞれの平均カウント:

cummean = np.cumsum(sample)/np.arange(1, n + 1)

は私がサイクルと理解することなく、中央値に対して同様の何かを行うことができますか?

答えて

2

ここでは、要素を行に沿って複製して2Dの配列を与えるアプローチがあります。次に、上の三角形の領域を大きな数字で塗りつぶして、後で各行に沿って配列を並べ替えると、基本的に対角要素まですべての要素をソートし、累積ウィンドウをシミュレートします。次に、真ん中または中間の2つの要素を選択するの定義に従って、最初の位置:(0,0)、次に2番目の行の要素を取得します。平均値は(1,0) & (1,1)です。第3の行について:(2,1)、第4の行について:(3,1) & (3,2)などの平均。したがって、ソートされた配列からそれらの要素を抽出し、したがって中央値を持ちます。

したがって、実装は次のようになり -

def cummedian_sorted(a): 
    n = a.size 
    maxn = a.max()+1 
    a_tiled_sorted = np.tile(a,n).reshape(-1,n) 
    mask = np.triu(np.ones((n,n),dtype=bool),1) 

    a_tiled_sorted[mask] = maxn 
    a_tiled_sorted.sort(1) 

    all_rows = a_tiled_sorted[np.arange(n), np.arange(n)//2].astype(float) 
    idx = np.arange(1,n,2) 
    even_rows = a_tiled_sorted[idx, np.arange(1,1+(n//2))] 
    all_rows[idx] += even_rows 
    all_rows[1::2] /= 2.0 
    return all_rows 

ランタイム試験

アプローチ -

# Loopy solution from @Uriel's soln 
def cummedian_loopy(arr): 
    return [median(a[:i]) for i in range(1,len(a)+1)] 

# Nan-fill based solution from @Nickil Maveli's soln 
def cummedian_nanfill(arr): 
    a = np.tril(arr).astype(float) 
    a[np.triu_indices(a.shape[0], k=1)] = np.nan 
    return np.nanmedian(a, axis=1) 

タイミング -

セット#1:

In [43]: a = np.random.randint(0,100,(100)) 

In [44]: print np.allclose(cummedian_loopy(a), cummedian_sorted(a)) 
    ...: print np.allclose(cummedian_loopy(a), cummedian_nanfill(a)) 
    ...: 
True 
True 

In [45]: %timeit cummedian_loopy(a) 
    ...: %timeit cummedian_nanfill(a) 
    ...: %timeit cummedian_sorted(a) 
    ...: 
1000 loops, best of 3: 856 µs per loop 
1000 loops, best of 3: 778 µs per loop 
10000 loops, best of 3: 200 µs per loop 

セット#2:

In [46]: a = np.random.randint(0,100,(1000)) 

In [47]: print np.allclose(cummedian_loopy(a), cummedian_sorted(a)) 
    ...: print np.allclose(cummedian_loopy(a), cummedian_nanfill(a)) 
    ...: 
True 
True 

In [48]: %timeit cummedian_loopy(a) 
    ...: %timeit cummedian_nanfill(a) 
    ...: %timeit cummedian_sorted(a) 
    ...: 
10 loops, best of 3: 118 ms per loop 
10 loops, best of 3: 47.6 ms per loop 
100 loops, best of 3: 18.8 ms per loop 

セット#3:

In [49]: a = np.random.randint(0,100,(5000)) 

In [50]: print np.allclose(cummedian_loopy(a), cummedian_sorted(a)) 
    ...: print np.allclose(cummedian_loopy(a), cummedian_nanfill(a)) 

True 
True 

In [54]: %timeit cummedian_loopy(a) 
    ...: %timeit cummedian_nanfill(a) 
    ...: %timeit cummedian_sorted(a) 
    ...: 
1 loops, best of 3: 3.36 s per loop 
1 loops, best of 3: 583 ms per loop 
1 loops, best of 3: 521 ms per loop 
+0

私は 'loopy'の方が速くなっていますが、統計情報ではなく' np.median'を使っていました。 – hpaulj

2

使用statistics.medianとcummulativeリスト内包(奇数インデックスが偶数長さリストの中央値が含まれていることに注意 - 中央値は、二つのメジアン要素の平均であり、それは通常、小数ではなく整数で生じる場合):

>>> from statistics import median 
>>> arr = [1, 3, 4, 2, 5, 3, 6] 
>>> cum_median = [median(arr[:i+1]) for i in range(len(arr)-1)] 
>>> cum_median 
[1, 2.0, 3, 2.5, 3, 3.0] 
+0

オープンサイクルと理解なし –

+0

'np.median'は' statistics.median'よりもはるかに高速です。 – hpaulj

1

は後半エントリの余地はありますか?

def cummedian_partition(a): 
    n = len(a) 
    assert n%4 == 0 # for simplicity 
    mn = a.min() - 1 
    mx = a.max() + 1 
    h = n//2 
    N = n + h//2 
    work = np.empty((h, N), a.dtype) 
    work[:, :n] = a 
    work[:, n] = 2*mn - a[0] 
    i, j = np.tril_indices(h, -1) 
    work[i, n-1-j] = (2*mn - a[1:h+1])[j] 
    k, l = np.ogrid[:h, :h//2 - 1] 
    work[:, n+1:] = np.where(k > 2*l+1, mx, 2 * mn - mx) 
    out = np.partition(work, (N-n//2-1, N-n//2, h//2-1, h//2), axis=-1) 
    out = np.r_[2*mn-out[:, h//2: h//2-2:-1], out[::-1, N-n//2-1:N-n//2+1]] 
    out[::2, 0] = out[::2, 1] 
    return np.mean(out, axis=-1) 

アルゴリズムは、線形の複雑さを持つパーティションを使用します。 np.partitionは1行ごとの分割点をサポートしていないため、体操が必要です。複雑さと必要なメモリの合計は二次的です。現在のベストに比べ

タイミング:

for j in (100, 1000, 5000): 
    a = np.random.randint(0, 100, (j,)) 
    print('size', j) 
    print('correct', np.allclose(cummedian_partition(a), cummedian_sorted(a))) 
    print('Divakar', timeit(lambda: cummedian_sorted(a), number=10)) 
    print('PP', timeit(lambda: cummedian_partition(a), number=10)) 

# size 100 
# correct True 
# Divakar 0.0022412699763663113 
# PP 0.002393342030700296 
# size 1000 
# correct True 
# Divakar 0.20881508802995086 
# PP 0.10222102201078087 
# size 5000 
# correct True 
# Divakar 6.158387024013791 
# PP 3.437395485001616 
2

Pythonはあなたが反復可能なためランニング「最小」を保つことができますheapqモジュールを持っていることを知って、私はheapqmedianで検索をした、と様々なアイテムを見つけましたsteaming medium。これ:

http://www.ardendertat.com/2011/11/03/programming-interview-questions-13-median-of-integer-stream/

は、値の下半分を有するもの、上半分と他の2つのheapqを維持class streamMedianを有しています。中央値は、1の「トップ」または両方の値の平均のいずれかです。クラスはinsertメソッドとgetMedianメソッドを持っています。ほとんどの作業はinsertです。

私はIpythonセッションにあることをコピーされ、定義された:

def cummedian_stream(b): 
    S=streamMedian() 
    ret = [] 
    for item in b: 
     S.insert(item) 
     ret.append(S.getMedian()) 
    return np.array(ret) 

テスト:

In [155]: a = np.random.randint(0,100,(5000)) 
In [156]: amed = cummedian_stream(a) 
In [157]: np.allclose(cummedian_sorted(a), amed) 
Out[157]: True 
In [158]: timeit cummedian_sorted(a) 
1 loop, best of 3: 781 ms per loop 
In [159]: timeit cummedian_stream(a) 
10 loops, best of 3: 39.6 ms per loop 

heapqストリームアプローチが道高速です。


@Urielが与えたリストの理解は比較的遅いです。私はそれが@Divakar'sソートソリューションよりも高速であるstatistics.medianためnp.medianに置き換えた場合でも:

def fastloop(a): 
    return np.array([np.median(a[:i+1]) for i in range(len(a))]) 

In [161]: timeit fastloop(a) 
1 loop, best of 3: 360 ms per loop 

そして@Paul Panzer'sパーティションのアプローチも良いですが、それでもストリーミングクラスに比べて遅くなります。

In [165]: timeit cummedian_partition(a) 
1 loop, best of 3: 391 ms per loop 

(必要に応じてstreamMedianクラスをコピーすることができます)。

+0

「numpy」でさえ、 'O(n^2)'(またはそれより悪い)のアルゴリズムを、それは 'O(n log n)'よりも速くすることができます。 –

関連する問題