这是
rvs从https://github.com/scipy/scipy/pull/5622/files中提取的方法，只需进行最小的更改-
足以作为独立的numpy函数运行。
import numpy as npdef rvs(dim=3): random_state = np.random 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()) # Householder transformation Hx = (np.eye(dim-n+1) - 2.*np.outer(x, x)/(x*x).sum()) mat = np.eye(dim) mat[n-1:, n-1:] = Hx H = np.dot(H, mat) # Fix the last sign such that the determinant is 1 D[-1] = (-1)**(1-(dim % 2))*D.prod() # Equivalent to np.dot(np.diag(D), H) but faster, apparently H = (D*H.T).T return H
它符合沃伦的测试，https://stackoverflow.com/a/38426572/901925
