function [sh2 ph2] = sbm_sample_h2(theta, sh1)

    n = size(sh1, 1);
            
    % only contributions from layer 1 hiddens and layer 2 biases
    raw = sh1 * theta.W2' + repmat_fast(theta.b2', n);

    % draw
    ph2 = sigmoid(raw/theta.T);
    sh2 = draw_bernoulli(ph2);
    
end