numpy配列のモードを得る最も速い方法

Question

私はhdf5ファイルから読み込んだNumPy配列のモードを見つけなければなりません。 NumPy配列は1dであり、浮動小数点値を含んでいます。

my_array=f1[ds_name].value  
mod_value=scipy.stats.mode(my_array) 

私の配列は1dで約1Mの値を含んでいます。私のスクリプトがモード値を返すのに約15分かかります。これをもっと速くする方法はありますか？

もう1つの質問は、モードが動作しているときにscipy.stats.median(my_array)が機能しない理由ですか？

AttributeError: module 'scipy.stats' has no attribute 'median'

Answer 1

scipy.stats.modeの実装では、多次元配列でaxis引数を処理するためのPythonのループを有します。

def mode1(x): 
    values, counts = np.unique(x, return_counts=True) 
    m = counts.argmax() 
    return values[m], counts[m] 

次の例は、1次元配列の場合の例です。まず、scipy.stats.modeとmode1は、同じ結果を与えることを長さ1000000

In [40]: x = np.random.randint(0, 1000, size=(2, 1000000)).sum(axis=0) 

In [41]: x.shape 
Out[41]: (1000000,) 

チェックして整数の配列を作ります。

In [42]: from scipy.stats import mode 

In [43]: mode(x) 
Out[43]: ModeResult(mode=array([1009]), count=array([1066])) 

In [44]: mode1(x) 
Out[44]: (1009, 1066) 

パフォーマンスを確認します。 mode(x)ため

In [45]: %timeit mode(x) 
2.91 s ± 18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) 

In [46]: %timeit mode1(x) 
39.6 ms ± 83.8 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) 

2.91秒とmode1(x)のための唯一の39.6ミリ秒。それはそれをソートするように、この変異するには/入力配列を変更すること

def mode1d(ar_sorted): 
    ar_sorted.sort() 
    idx = np.flatnonzero(ar_sorted[1:] != ar_sorted[:-1]) 
    count = np.empty(idx.size+1,dtype=int) 
    count[1:-1] = idx[1:] - idx[:-1] 
    count[0] = idx[0] + 1 
    count[-1] = ar_sorted.size - idx[-1] - 1 
    argmax_idx = count.argmax() 

    if argmax_idx==len(idx): 
     modeval = ar_sorted[-1] 
    else: 
     modeval = ar_sorted[idx[argmax_idx]] 
    modecount = count[argmax_idx] 
    return modeval, modecount 

注 -

Answer 2

ここでのソートに基づいて一つの手法です。したがって、入力配列をに変更せずに保存するか、ソートされている入力配列が気になる場合は、コピーを渡してください。 1M要素に

サンプル実行 - コピーして

In [65]: x = np.random.randint(0, 1000, size=(1000000)).astype(float) 

In [66]: from scipy.stats import mode 

In [67]: mode(x) 
Out[67]: ModeResult(mode=array([ 295.]), count=array([1098])) 

In [68]: mode1d(x) 
Out[68]: (295.0, 1098) 

ランタイムテスト

In [75]: x = np.random.randint(0, 1000, size=(1000000)).astype(float) 

# Scipy's mode 
In [76]: %timeit mode(x) 
1 loop, best of 3: 1.64 s per loop 

# @Warren Weckesser's soln 
In [77]: %timeit mode1(x) 
10 loops, best of 3: 52.7 ms per loop 

# Proposed in this post 
In [78]: %timeit mode1d(x) 
100 loops, best of 3: 12.8 ms per loop 

、mode1dのタイミングはmode1に匹敵するだろう。

Answer 3

上記の返信の2つの関数mode1とmode1dをスクリプトに追加し、scipy.stats.modeと比較しようとしました。

dir_name="C:/Users/test_mode" 
file_name="myfile2.h5" 
ds_name="myds" 
f_in=os.path.join(dir_name,file_name) 

def mode1(x): 
    values, counts = np.unique(x, return_counts=True) 
    m = counts.argmax() 
    return values[m], counts[m] 

def mode1d(ar_sorted): 
    ar_sorted.sort() 
    idx = np.flatnonzero(ar_sorted[1:] != ar_sorted[:-1]) 
    count = np.empty(idx.size+1,dtype=int) 
    count[1:-1] = idx[1:] - idx[:-1] 
    count[0] = idx[0] + 1 
    count[-1] = ar_sorted.size - idx[-1] - 1 
    argmax_idx = count.argmax() 

    if argmax_idx==len(idx): 
     modeval = ar_sorted[-1] 
    else: 
     modeval = ar_sorted[idx[argmax_idx]] 
    modecount = count[argmax_idx] 
    return modeval, modecount 


startTime=time.time() 
with h5py.File(f_in, "a") as f1: 

     myds=f1[ds_name].value 
     time1=time.time() 
     file_read_time=time1-startTime 
     print(str(file_read_time)+"\t"+"s"+"\t"+str((file_read_time)/60)+"\t"+"min") 

     print("mode_scipy=") 
     mode_scipy=scipy.stats.mode(myds) 
     print(mode_scipy) 
     time2=time.time() 
     mode_scipy_time=time2-time1 
     print(str(mode_scipy_time)+"\t"+"s"+"\t"+str((mode_scipy_time)/60)+"\t"+"min") 

     print("mode1=") 
     mode1=mode1(myds) 
     print(mode1) 
     time3=time.time() 
     mode1_time=time3-time2 
     print(str(mode1_time)+"\t"+"s"+"\t"+str((mode1_time)/60)+"\t"+"min") 

     print("mode1d=") 
     mode1d=mode1d(myds) 
     print(mode1d) 
     time4=time.time() 
     mode1d_time=time4-time3 
     print(str(mode1d_time)+"\t"+"s"+"\t"+str((mode1d_time)/60)+"\t"+"min") 

1Mの周囲のnumpyのアレイのためのスクリプトを実行した結果である：

mode_scipy = ModeResult（モード=配列（[1.11903353e-06]、DTYPE =のfloat32）、カウント=配列（[304909]）） 938.8368742465973 S
15.647281237443288分

MODE1 =（1。1190335e-06、304909）

0.06500649452209473 S
0.0010834415753682455分

mode1d =（1.1190335e-06、304909）

0.06200599670410156 S
0.0010334332784016928分