Bernoulli random number generator

Question

I cannot understand how Bernoulli Random Number generator used in numpy is calculated and would like some explanation on it. For example:

np.random.binomial(size=3, n=1, p= 0.5)

Results:
[1 0 0]

n = number of trails

p = probability of occurrence

size = number of experiments

With how do I determine the generated numbers/results of "0" or "1"?

=================================Update==================================

I created a Restricted Boltzmann Machine which always presents the same results despite being "random" on multiple code executions. The randomize is seeded using np.random.seed(10)

import numpy as np

np.random.seed(10)

def sigmoid(u):
    return 1/(1+np.exp(-u))

def gibbs_vhv(W, hbias, vbias, x):
    f_s = sigmoid(np.dot(x, W) + hbias)
    h_sample = np.random.binomial(size=f_s.shape, n=1, p=f_s)

    f_u = sigmoid(np.dot(h_sample, W.transpose())+vbias)
    v_sample = np.random.binomial(size=f_u.shape, n=1, p=f_u)
    return [f_s, h_sample, f_u, v_sample]

def reconstruction_error(f_u, x):
    cross_entropy = -np.mean(
        np.sum(
            x * np.log(sigmoid(f_u)) + (1 - x) * np.log(1 - sigmoid(f_u)),
            axis=1))
    return cross_entropy


X = np.array([[1, 0, 0, 0]])

#Weight to hidden
W = np.array([[-3.85, 10.14, 1.16],
              [6.69, 2.84, -7.73],
              [1.37, 10.76, -3.98],
              [-6.18, -5.89, 8.29]])

hbias = np.array([1.04, -4.48, 2.50]) #<= 3 bias for 3 neuron in hidden
vbias = np.array([-6.33, -1.68, -1.25, 3.45]) #<= 4 bias for 4 neuron in input


k = 2
v_sample = X
for i in range(k):
    [f_s, h_sample, f_u, v_sample] = gibbs_vhv(W, hbias, vbias, v_sample)
    start = v_sample
    if i < 2:
        print('f_s:', f_s)
        print('h_sample:', h_sample)
        print('f_u:', f_u)
        print('v_sample:', v_sample)
    print(v_sample)
    print('iter:', i, ' h:', h_sample, ' x:', v_sample, ' entropy:%.3f'%reconstruction_error(f_u, v_sample))

Results:

[[1 0 0 0]]
f_s: [[ 0.05678618  0.99652957  0.97491304]]
h_sample: [[0 1 1]]
f_u: [[ 0.99310473  0.00139984  0.99604968  0.99712837]]
v_sample: [[1 0 1 1]]
[[1 0 1 1]]
iter: 0  h: [[0 1 1]]  x: [[1 0 1 1]]  entropy:1.637

f_s: [[  4.90301318e-04   9.99973278e-01   9.99654440e-01]]
h_sample: [[0 1 1]]
f_u: [[ 0.99310473  0.00139984  0.99604968  0.99712837]]
v_sample: [[1 0 1 1]]
[[1 0 1 1]]
iter: 1  h: [[0 1 1]]  x: [[1 0 1 1]]  entropy:1.637

Answer 1

I am asking on how the algorithm works to produce the numbers. – WhiteSolstice 35 mins ago

Non-technical explanation

If you pass n=1 to the Binomial distribution it is equivalent to the Bernoulli distribution. In this case the function could be thought of simulating coin flips. size=3 tells it to flip the coin three times and p=0.5 makes it a fair coin with equal probabilitiy of head (1) or tail (0).

The result of [1 0 0] means the coin came down once with head and twice with tail facing up. This is random, so running it again would result in a different sequence like [1 1 0], [0 1 0], or maybe even [1 1 1]. Although you cannot get the same number of 1s and 0s in three runs, on average you would get the same number.

Technical explanation

Numpy implements random number generation in C. The source code for the Binomial distribution can be found here. Actually two different algorithms are implemented.

• If n * p <= 30 it uses inverse transform sampling.
• If n * p > 30 the BTPE algorithm of (Kachitvichyanukul and Schmeiser 1988) is used. (The publication is not freely available.)

I think both methods, but certainly the inverse transform sampling, depend on a random number generator to produce uniformly distributed random numbers. Numpy internally uses a Mersenne Twister pseudo random number generator. The uniform random numbers are then transformed into the desired distribution.

Answer 2

A Binomially distributed random variable has two parameters n and p, and can be thought of as the distribution of the number of heads obtained when flipping a biased coin n times, where the probability of getting a head at each flip is p. (More formally it is a sum of independent Bernoulli random variables with parameter p).

For instance, if n=10 and p=0.5, one could simulate a draw from Bin(10, 0.5) by flipping a fair coin 10 times and summing the number of times that the coin lands heads.

In addition to the n and p parameters described above, np.random.binomial has an additional size parameter. If size=1, np.random.binomial computes a single draw from the Binomial distribution. If size=k for some integer k, k independent draws from the same Binomial distribution will be computed. size can also be an array of indices, in which case a whole np.array with the given size will be filled with independent draws from the Binomial distribution.

Note that the Binomial distribution is a generalisation of the Bernoulli distribution - in the case that n=1, Bin(n,p) has the same distribution as Ber(p).

For more information about the binomial distribution see: https://en.wikipedia.org/wiki/Binomial_distribution