code: def rmv_normal_chol(mu, sig, size=1):
      
      mu_size = np.shape(mu)
  
      if size == 1:
          out = np.random.normal(size=mu_size)
          try:
              flib.dtrmm_wrap(sig, out, 'L', 'N', 'L', 1.)
          except:
              out = np.dot(sig, out)
          out += mu
          return out
      else:
          if not hasattr(size, '__iter__'):
              size = (size,)
          tot_size = np.prod(size)
          out = np.random.normal(size=(tot_size,) + mu_size)
          for i in xrange(tot_size):
              try:
                  flib.dtrmm_wrap(sig, out[i, :], 'L', 'N', 'L', 1.)
              except:
                  out[i, :] = np.dot(sig, out[i,:])
              out[i, :] += mu
          return out.reshape(size + mu_size)
docs: Random multivariate normal variates.