statistics.py 文件源码

def sample_(self, shape, a):
        if np.isscalar(a) or a.size == 1:
            return self.sample_scalar(shape, a)

        #print shape, a.shape
        AMAX = 30
        k = 1
        K = np.full(shape, k)

        # for large a values, just use non-truncated poisson
        large = broadcast_mask(a > AMAX)
        small = broadcast_mask(a <= AMAX)
        K[large] = np.random.poisson(a[large], K[large].shape)
        Ksmall = K[small]
        a = a[small]

        # actual algorithm begins here
        s = a/np.expm1(a)
        S = s
        U = np.random.random(Ksmall.shape)
        new = S < U
        while np.any(new):
            k += 1
            Ksmall[new] = k
            s = s*a/float(k)
            S = S + s
            new = S < U
        K[small] = Ksmall
        return K