Nothing
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
################################################################################
# FUNCTION: DESCRIPTION:
# isPositiveDefinite Checks if the matrix X is positive definite
# makePositiveDefinite Forces the matrix x to be positive definite
################################################################################
isPositiveDefinite <-
function(x)
{
# A function implemented by Diethelm Wuertz
# Description:
# Checks if the matrix x is positive definite
# Arguments:
# x - a symmetric matrix or any other rectangular object
# describing a covariance matrix which can be converted into
# a matrix by the function 'as.matrix'.
# FUNCTION:
# Transform:
x = as.matrix(x)
# Check if matrix is positive definite:
ans = .is.positive.definite(m = x)
# Return Value:
ans
}
# ------------------------------------------------------------------------------
.is.positive.definite <-
function(m, tol, method = c("eigen", "chol"))
{
# Author:
# Copyright 2003-05 Korbinian Strimmer
# Rank, condition, and positive definiteness of a matrix
# GNU General Public License, Version 2
method = match.arg(method)
if (!is.matrix(m)) {
m = as.matrix(m)
}
if (method == "eigen") {
eval = eigen(m, only.values = TRUE)$values
if( missing(tol) ) {
tol = max(dim(m))*max(abs(eval))*.Machine$double.eps
}
if (sum(eval > tol) == length(eval)) {
return(TRUE)
} else {
return(FALSE)
}
} else if (method == "chol") {
val = try(chol(m), silent = TRUE)
if (inherits(val, "try-error")) {
return(FALSE)
} else {
return(TRUE)
}
}
}
# ------------------------------------------------------------------------------
makePositiveDefinite <-
function(x)
{
# A function implemented by Diethelm Wuertz
# Description:
# Forces the matrix x to be positive definite
# Arguments:
# x - a symmetric matrix or any other rectangular object
# describing a covariance matrix which can be converted into
# a matrix by the function 'as.matrix'.
# FUNCTION:
# Make Positive Definite:
ans = .make.positive.definite(m = x)
# Return Value:
ans
}
# ------------------------------------------------------------------------------
.make.positive.definite <-
function(m, tol)
{
# Author:
# Copyright 2003-05 Korbinian Strimmer
# Rank, condition, and positive definiteness of a matrix
# GNU General Public License, Version 2
# Method by Higham 1988
if (!is.matrix(m)) {
m = as.matrix(m)
}
d = dim(m)[1]
if ( dim(m)[2] != d ) {
stop("Input matrix is not square!")
}
es = eigen(m)
esv = es$values
if (missing(tol)) {
tol = d*max(abs(esv))*.Machine$double.eps
}
delta = 2*tol
# factor two is just to make sure the resulting
# matrix passes all numerical tests of positive definiteness
tau = pmax(0, delta - esv)
dm = es$vectors %*% diag(tau, d) %*% t(es$vectors)
# print(max(DA))
# print(esv[1]/delta)
return( m + dm )
}
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.