Compute orthonormal basis for the column span of X.
Rank is determined by zeroing all singular values, u, less than or equal to tol*u.max().
Adapted from spm_reml.m
ReML estimation of covariance components from sigma using design matrix.
sigma – m-by-m covariance matrix components – q-by-m-by-m array of variance components
mean of sigma is modeled as a some over components[i]
n – degrees of freedom of sigma penalty_cov – quadratic penalty to be applied in Fisher algorithm.
If the value is a float, f, the penalty is f * identity(m). If the value is a 1d array, this is the diagonal of the penalty.
C – estimated mean of sigma h – array of length q representing coefficients
of variance components
cov_h – estimated covariance matrix of h