#はじめに
Qiita初投稿です。よろしくお願い致します。
いきなり本題だが、
二項分布を用いたモデリングでの尤度計算にて、
組み合わせの計算 $nCr$ がたくさん登場するらしい。
ところで、0から8までしか値を取らないデータが1000件あって、
それぞれで
$n=8$ の組み合わせ計算 $nCr$ を行うのは、
無駄なのでは？
同じ計算を100回以上行わせてしまうのでは？
と考えた。
そこで、Pythonの辞書を扱い、
「あらかじめ全体の数$N$までの
全ての整数における組み合わせを計算しておくのはどうだろう？」
と思い、何パターンか時間計測をしてみた。
また初めてのQiitaへの投稿の練習もさせていただきたい。
#組み合わせ計算はscipy.specialを使用
@derodero24 様がまとめてくださった、
【Python】組み合わせ(nCr) 計算の高速化
(https://qiita.com/derodero24/items/91b6468e66923a87f39f)
を参考に、
scipy.special.comb()を使用する。
from scipy.special import comb
a=comb(n,r)
また、二項分布の乱数発生は
scipy.stats.binom()を使用。
発生した乱数 $r_i$ それぞれについて、
$nCr_i$ を求め、
辞書未使用と使用でどれだけ計算時間がかかるか
比較していく。
#未使用パターン
（マイブームなので）今回は、numpyのfrompyfuncを使用し、
scipy.special.comb()をグローバル関数化する。
import numpy as np
from scipy.stats import binom
from scipy.special import comb
def comb_frompyfunc(N,d_size,iter_count):
array_list=[]
comb_np=np.frompyfunc(comb,2,1)
for _ in range(iter_count):
d=binom.rvs(N,p,size=d_size)
comb_array=comb_np(N,d)
array_list.append(comb_array)
return array_list
(numpy.vectorize(comb)(N,d)などでも計算できるが、今回は使用しない。また、事前に比較したところ、frompyfuncを使ったほうが、どのパターンでも早かった)
#辞書使用パターン
二項分布による乱数を発生させる前に、
辞書comb_dictに、0からNまでの組み合わせの数を代入し、
np.vectorize(comb_dict.get)(d)
で乱数それぞれの組み合わせの数に変換する。
import numpy as np
from scipy.stats import binom
from scipy.special import comb
def comb_fromdict(N,d_size,iter_count):
array_list=[]
comb_dict={}
for i in range(0,N+1):
comb_dict[i]=comb(N,i)
for _ in range(iter_count):
d=binom.rvs(N,p,size=d_size)
comb_array=np.vectorize(comb_dict.get)(d)
array_list.append(comb_array)
return array_list
辞書の作成方法について、forループを使わずに行う方法を知っていればよかったが、上手く行かなかったり、逆に遅くなったりしたので、今回は妥協。
#計測パターン１ 組み合わせ数 ＜ データサイズ
組み合わせ $N=10$
データサイズ $d\_size=100$
N=10
d_size=100
iter_count=1
%timeit comb_frompyfunc(N,d_size,iter_count) #未使用
%timeit comb_fromdict(N,d_size,iter_count) #辞書使用
"""
出力
1000 loops, best of 3: 751 µs per loop
1000 loops, best of 3: 196 µs per loop
"""
辞書使用のほうが3.8倍早くなった。
#計測パターン2 組み合わせ数 = データサイズ
組み合わせ $N=10$
データサイズ $d\_size=10$
N=10
d_size=10
iter_count=1
%timeit comb_frompyfunc(N,d_size,iter_count) #未使用
%timeit comb_fromdict(N,d_size,iter_count) #辞書使用
"""
出力
10000 loops, best of 3: 123 µs per loop
10000 loops, best of 3: 158 µs per loop
"""
おおかた予想通りだが、
組み合わせ数とデータサイズが等しい場合は、
辞書を使用したほうが21%ほど遅い。
次は組み合わせ数もデータサイズも1000に変更してみた
組み合わせ $N=1000$
データサイズ $d\_size=1000$
N=1000
d_size=1000
iter_count=1
%timeit comb_frompyfunc(N,d_size,iter_count) #未使用
%timeit comb_fromdict(N,d_size,iter_count) #辞書使用
"""
出力
100 loops, best of 3: 6.26 ms per loop
100 loops, best of 3: 6.52 ms per loop
"""
今回も辞書型のほうが遅かったが、
約4%の差に減少している。
#計測パターン3 組み合わせ数 <<< データサイズ
組み合わせ $N=10$
データサイズ $d\_size=10000$
N=10
d_size=10000
iter_count=1
%timeit comb_frompyfunc(N,d_size,iter_count) #未使用
%timeit comb_fromdict(N,d_size,iter_count) #辞書使用
"""
出力
10 loops, best of 3: 59 ms per loop
1000 loops, best of 3: 1.72 ms per loop
"""
辞書使用のほうが34.3倍早くなった。
計測パターン１と比較して、
データサイズが1000倍にはなったが、
時間計測の差は10倍程度しか広がらず。
#まとめ
組み合わせ数含めた計算は、データの量が大きくなればなるほど
pythonの辞書を使用したほうが、早くなることがわかった。
辞書の作成にforループを使う以外で適した方法があれば更に早くなるような気はする。
#numpyによる配列の参照-2019/08/15追記
投稿後、@konandoiruasa様より、
コメントにより、numpy.arrayによる方法で、断然早くなることを
ご教授いただきました。
def comb_fromtable(N,d_size,iter_count):
comb_table = np.array([comb(N, i) for i in range(N+1)])
d = binom.rvs(N, p, size=d_size*iter_count)
return comb_table[d]
ありがとうございます！
numpy.arrayで組み合わせ配列を作ってしまい、
comb_table[d]と指定してしまえば、
私のやりたかったことが一気に実現してしまいました。
また、size = d_size * iter_count と指定してしまえば、
forループやappendをする必要もありませんでした。
これでまた2つ賢くなりました。
Qiitaに記事を投稿してみてよかった瞬間でございました。