numpyのループでこれをベクトル化する

Question

私はnumpyでPythonで行列z（以下に定義）を計算しようとしています。

ここで（forループの1を使用して）私の現在のソリューション

z = np.zeros((n, k)) 
for i in range(n): 
    v = pi * (1/math.factorial(x[i])) * np.exp(-1 * lamb) * (lamb ** x[i]) 
    numerator = np.sum(v) 
    c = v/numerator 
    z[i, :] = c 
return z 

が、それは完全にこの計算をベクトル化することは可能ですね。私は数千回の反復でこの計算を行う必要があり、numpyでの行列演算はforループの場合よりもはるかに高速です。

Answer 1

ここにはEのベクトル化バージョンがあります。それはNumPy broadcastingとアレイベースの算術for-loopとスカラー算術を置き換える：

def alt_E(x): 
    x = x[:, None] 
    z = pi * (np.exp(-lamb) * (lamb**x))/special.factorial(x) 
    denom = z.sum(axis=1)[:, None] 
    z /= denom 
    return z 

---

私はx、lamb、pi、nとkの典型的なサイズのための感覚を得るためにem.pyを走りました。このサイズのデータ​​では、 alt_EはE比べて約120xの速いです：

In [32]: %timeit E(x) 
100 loops, best of 3: 11.5 ms per loop 

In [33]: %timeit alt_E(x) 
10000 loops, best of 3: 94.7 µs per loop 

In [34]: 11500/94.7 
Out[34]: 121.43611404435057 

---

これは私がベンチマークに使用する設定です：ところで

import math 
import numpy as np 
import scipy.special as special 

def alt_E(x): 
    x = x[:, None] 
    z = pi * (np.exp(-lamb) * (lamb**x))/special.factorial(x) 
    denom = z.sum(axis=1)[:, None] 
    z /= denom 
    return z 


def E(x): 
    z = np.zeros((n, k)) 
    for i in range(n): 
     v = pi * (1/math.factorial(x[i])) * \ 
      np.exp(-1 * lamb) * (lamb ** x[i]) 
     numerator = np.sum(v) 
     c = v/numerator 
     z[i, :] = c 
    return z 

n = 576 
k = 2 

x = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 5]) 

lamb = np.array([ 0.84835141, 1.04025989]) 
pi = np.array([ 0.5806958, 0.4193042]) 

assert np.allclose(alt_E(x), E(x)) 

---

、Escipy.stats.poissonを使用して計算することもできます。

import scipy.stats as stats 
pois = stats.poisson(mu=lamb) 

def alt_E2(x): 
    z = pi * pois.pmf(x[:,None]) 
    denom = z.sum(axis=1)[:, None] 
    z /= denom 
    return z 

が、これは、少なくともこの長さの配列のために、より高速であることが判明していません。

In [33]: %timeit alt_E(x) 
10000 loops, best of 3: 94.7 µs per loop 

In [102]: %timeit alt_E2(x) 
1000 loops, best of 3: 278 µs per loop 

をxより大きいため、alt_E2が高速です：

In [104]: x = np.random.random(10000) 

In [106]: %timeit alt_E(x) 
100 loops, best of 3: 2.18 ms per loop 

In [105]: %timeit alt_E2(x) 
1000 loops, best of 3: 643 µs per loop