#はじめに

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()を使用する。

sample.py

from scipy.special import comb

a=comb(n,r)

また、二項分布の乱数発生は
scipy.stats.binom()を使用。

発生した乱数 $r_i$ それぞれについて、
$nCr_i$ を求め、
辞書未使用と使用でどれだけ計算時間がかかるか
比較していく。

#未使用パターン

（マイブームなので）今回は、numpyのfrompyfuncを使用し、
scipy.special.comb()をグローバル関数化する。

comb_frompyfunc(N,d_size,iter_count)

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)
で乱数それぞれの組み合わせの数に変換する。

comb_fromdict(N,d_size,iter_count)

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$

pattern1.py

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$

pattern2_1.py

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$

pattern2_2.py

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$

pattern3.py

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に記事を投稿してみてよかった瞬間でございました。