'Draw random samples from SO(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)
    D = np.ones((dim,))
    for n in range(1, dim):
        x = random_state.normal(size=(((dim - n) + 1),))
        D[(n - 1)] = np.sign(x[0])
        x[0] -= (D[(n - 1)] * np.sqrt((x * x).sum()))
        Hx = (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)
    D[(-1)] = (((-1) ** (1 - (dim % 2))) * D.prod())
    H = (D * H.T).T
    return H
