tr = approximatetraceof_inverse(A)

bentB = bending(B [, W], gamma=nothing, tol=eps, verbose=false)
bentB = bending(B, gamma=nothing, tol=eps, verbose=false)

Apply a traditional "bending" method to a symmetric matrix B (between-group-variance) using another symmetric matrix W (within-group-variance), as described by Hayes and Hill (1981). The resulting matrix is bentB = (1-gamma)*B + gamma*v*W, where v is the average eigenvalues of B and gamma is a constant between 0 and 1, which can be provided by the user or generated by this function. When W is omitted, the identity matrix will be used instead.

If the user supplies gamma, the function will simply apply the bending formula and return bentB without checking whether bentB is positive defnite or not. Otherwise (gamma=nothing), this function generates gamma=1/1000, applies bending, and returns bentB if it is positive definite as minimum(eigen(bentB).values)>tol (the default tol is the machine epsilon). If not, the function tries gamma=2/1000, gamma=3/1000, ..., and gamma=1000/1000 until bentB becomes positive definite.

With verbose=true, this function prints v (average eigenvalue of B) and gamma that has been used for bending (default: verbose=false).

bentV = bending2(V [, W], corr=false, mineig=eps, tol=eps, maxiter=10000, force=false, verbose=false)

Apply an alternative "bending" method to a variance-covariance matrix V using a weighting matrix W, as described by Jorjani et al. (2003). When omitting W, J (matrix of ones) will be used, implying the same weights for all elements. Jorjani et al. (2003) suggested W as the reciprocal of number of individuals used in the estimation of covariances.

The function calculates the eigenvalues of V, and an eigenvalue is replaced with mineig if it is smaller than mineig. The covariance matrix is reconstructed with modified eigenvalues, then tests whether it is positive definite by checking the minimum eigenvalue larger than tol. The function repeats the process until the resulting matrix becomes positive definite, and the number of maximum iterations is defind by maxiter.

See some options below.

• With corr=true, this function assumes the input is a correlation matrix; otherwise, the input is expected to be a variance-covariance matrix.
• With force=true, this function iterates the formula maxiter times regardless of whether it has converged or not.
• With verbose=true, this function shows the number of iterations required when force=false.

 dv,ev = chol_symtrid(T::SymTridiagonal)

Calculate the Cholesky factor of a symmetric tridiagonal matrix. Note that the factor is a lower tridiagonal matrix, and the output is a couple of vectors. The factor is stored in dv for the diagonal matrix (D) and in ev for the off-diagonal elements (L).

 M = directsum(X...)

Returns the direct sum of supplied matrices.

disp(x)
disp(x...)
disp(io::IO,x)
disp(io::IO,x...)
disp(m::Array; format)
disp(io::IO, m::Array; format)

Write an object x (or objects x...) to io (or stdout if io is not given). The output is the same as println and Base.print_matrix, but the element type will not be shown. If an array is given, you can specify the C-style format with format.

A = [1 2; 3 4]
b = [9,10]

disp(A)
disp(b)
disp(A,b)

Z = get_design_matrix(x, nested=[]; maxlev=0, cov=false)
Z = get_design_matrix(x, y)
Z = get_design_matrix(X; maxlev=[], cov=[])

Returns a design matrix accoding to class variables or covariates in a vector x or matrix X. For a vector x, optional argument dested defines x as the nested covariate within a class. Option maxval specifies the number of levels in this effect. Omitting maxval, the function uses the maximum integer as the maximum level in x or X. Option cov specifies the input is covariate and passes it through the output. The matrix X doesn't support nested covariates.

If two vectors x and y are provided, the function returns the design matrix for the intercation between the two effects; the main effects are not included in the result. Also, in this case, the function does not accept any additional options.

julia> x = [1,2,1,1,2]
julia> X = get_design_matrix(x)
5×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
 1.0  0.0
 1.0  0.0
 0.0  1.0

X = get_interaction_matrix(X1, X2)

Returns a design matrix for the interaction between effects. The design matrices for the two main effects (X1 and X2) should be provided. Note that this function will not check if the supplied matrices are really the design matrices - it just creats all the combination of column products among the input matrices.

vare = hrv_class(ve,rangeset)

Generate a vector containing residual variances from a vector of heterogeneous residual variances ve and a set of ranges rangeset.

julia> vare = hrv_class([10.0,20.0],[3:5,6:9])
7-element Vector{Float64}:
10.0
10.0
10.0
20.0
20.0
20.0
20.0

dv,ev = ldlt_symtrid(T::SymTridiagonal)

Calculate the LDLT factor of a symmetric tridiagonal matrix. Note that the factor is a lower tridiagonal matrix, and the output is a couple of vectors. The factor is stored in dv for the diagonal matrix (D) and in ev for the off-diagonal elements (L).

L = legendre_polynomials(n; normalize=true)

Calculates an upper triangular matrix with the n-th order Legendre coefficients. By default, the coefficients will be normalized as normalize=true.

data = load_data_set(name)

Load a small data set for numerical testing.

v = normalized_legendre(x,n)

Calculates a single value of the n-th order Legendre polynomial of x. It calls Jacobi.legendre.

M = normalized_legendre_matrix(t,n; tmin=minimum(t),tmax=maximum(t))

Calculates a matrix with the n-th order Legendre coefficients on time points t, a vector or range. The earliest and latest (minimum and maximum) time points will be taken from the edges of t. You can provide these time points, tmin and tmax, as options.

# up to 3rd order coefficients for 5 to 12
julia> M = normalized_legendre_matrix(5:12,3)

# same expression
julia> M = normalized_legendre_matrix([5,6,7,8,9,10,11,12],3)

# as above with time points ranging 1 and 15
julia> M = normalized_legendre_matrix(5:12,3, tmin=1, tmax=15)

cov = r2v(corr, stdev)

Convert a correlation matrix cov and the vector of standard deviation stdev to the covariance matrix cov. This is a wrapper of StatsBase.cor2cov.

x = solve_pcg(A,b; maxiter=5000,eps=1e-15,reset=50,verbose=false)

Solve the symmetric positive definite system of equations: A as sparse left-hand side and b as a vector of right-hand side. The Jacobi preconditioner (diag(A)) is used. The maximum number of iterations is specified with maxiter. Convergence criterion eps is norm(b-A*x)/norm(b). The residual vector will be recalculated every reset rounds. The details will be shown with verbose=true.

julia> x = solve_pcg(A,b)

takahashi_ldlt_symtrid!(dv,ev)

Compute the sparse inverse (as known as selected inverse) of the LDLT factor of a symmetric trigiagonal matrix. The LDLT factor is stored in two vectors (dv and ev), which should be calculated with ldlt_symtrid.

A = unvech(v; n)

Convert a half-stored vector v to a square matrix A with order n. It is the inverse function of vech. When n is missing, the function will guess the order.

corr = v2r(cov)

Convert a covariance matrix cov to the correlation matrix corr. This is a wrapper of StatsBase.cov2cor.

v = vech(A)

Convert a square matrix to a half-stored vector. It was written by Steven G. Johnson; see https://discourse.julialang.org/t/half-vectorization/7399 for details.