'Draw random samples from O(N).
Parameters
dim : integer
Dimension of rotation space (N).
size : integer, optional
Number of samples to draw (default 1).
Returns
rvs : ndarray or scalar
Random size N-dimensional matrices, dimension (size, dim, dim)'

def rvs(self, dim, size=1, random_state=None):

    size = int(size)
    if (size > 1):
        return np.array([self.rvs(dim, size=1, random_state=random_state) for i in range(size)])
    dim = self._process_parameters(dim)
    random_state = self._get_random_state(random_state)
    H = np.eye(dim)
    for n in range(1, dim):
        x = random_state.normal(size=(((dim - n) + 1),))
        D = np.sign(x[0])
        x[0] += (D * np.sqrt((x * x).sum()))
        Hx = ((- D) * (np.eye(((dim - n) + 1)) - ((2.0 * np.outer(x, x)) / (x * x).sum())))
        mat = np.eye(dim)
        mat[(n - 1):, (n - 1):] = Hx
        H = np.dot(H, mat)
    return H
