code: def rskew_normal(mu, tau, alpha, size=()):
      
      size_ = size or (1,)
      len_ = np.prod(size_)
      return flib.rskewnorm(
          len_, mu, tau, alpha, np.random.normal(size=2 * len_)).reshape(size)
docs: Skew-normal random variates.