How hidden layer is made binary in Restricted Boltzmann Machine (RBM)?
As you correctly say, we calculate the probability of a hidden unit $h_j$
being one and then make it binary. That probability is given by
$$p(h_j=1) = \sigma\left(b_j + \sum_{i=1}^V w_{ij}v_i \right)$$
where $\sigma$ is the sigmoid function, $b_j$ is the bias of hidden unit $h_j$, $V$ is the number of visible units, $v_i$ is the (binary!) state of visible unit $i$, and $w_{ij}$ are the weights.

So, your MATLAB code for obtaining the probabilities hidden_probs is something like this (we write the sum implicitly by making a matrix multiplication):

hidden_probs = sigmoid(hidden_bias + data * weights)


Now, we have the probability $p(h_j=1)$ for each hidden unit $j \in [1,H]$. Now, this is only a probability. And we need a binary number, either 0 or 1. So the only thing we can do is pick a random sample from the probability distribution of $h_j$, which is a Bernoulli distribution.

As all hidden units are independent, we need to get one sample for each hidden unit independently. And also, in each training step, we need to draw new samples.

To draw these samples from the Bernoulli distribution, you can use the built-in functions of e.g. MATLAB (binornd) or Python (numpy.random.binomial). Note that these functions are to sample from a binomial distribution, but the Bernoulli distribution is just a special case of the binomial distribution with N=1. In MATLAB, that would be something like

hidden_states = binornd(1, hidden_probs)


which would create vector hidden_states which contains either 0 or 1, drawn randomly for each probability in hidden_probs.

As you probably have noticed, nobody does that! E.g. describes it in his Practical Guide to Training RBMs, as


  the hidden unit turns on if this probability is greater than a random number uniformly distributed between 0 and 1.


That is exactly what Hinton does in his RBM code: he gets a random number for each hidden unit using rand, i.e. randomly sampled from the uniform distribution between [0,1]. He then does the comparison:

hidden_states = hidden_probs &gt; rand(1, H)


This is equivalent to using binornd, but is probably faster. For example to generate a random number that is 1 with p=0.9, you get a random number from [0,1]. Now, in 90% of the cases, this random number is smaller than 0.9, and in 10% of the cases it is larger than 0.9. So to get a random number that is 1 with p=0.9, you can call 0.9 &gt; rand(1) - which is exactly what they do.

tl;dr: Pick a new random number from the range [0,1] for each hidden unit in each iteration. Compare it to your probability with hidden_probs &gt; rand(1,H) to make it binary.