sns.check.logdensity: Utility function for validating log-density

In sns: Stochastic Newton Sampler (SNS)

Utility function for validating log-density

Description

Utility function for validating log-density: 1) dimensional consistency of function argument, gradient and Hessian, 2) finiteness of function, gradient and Hessian, 3) closeness of analytical and numerical derivatives, and 4) negative definiteness of Hessian.

Usage

sns.check.logdensity(x, fghEval
  , numderiv.method = c("Richardson", "complex")
  , numderiv.args = list()
  , blocks = append(list(1:length(x)), as.list(1:length(x)))
  , dx = rep(1, length(x)), nevals = 100, negdef.tol = 1e-08, ...)
## S3 method for class 'sns.check.logdensity'
print(x, ...)

Arguments

x	For sns.check.logdensity, initial point, around which a random collection of points are generated to perform validation tests. For print.sns.check.logdensity, an object of class sns.check.logdensity, typically the output of sns.check.logdensity function.
fghEval	Log-density to be validated. A valid log-density can have one of 3 forms: 1) return log-density, but no gradient or Hessian, 2) return a list of f and g for log-density and its gradient vector, respectively, 3) return a list of f, g, and h for log-density, gradient vector, and Hessian matrix.
numderiv.method	Method used for numeric differentiation. This is passed to the grad and hessian functions in numDeriv package. See the package documentation for details.
numderiv.args	Arguments to the numeric differentiation method chosen in numderiv.method, passed to grad and hessian functions in numDeriv. See package documentation for details.
blocks	A list of state space subsets (identified by their positional indexes), for which negative-definiteness of Hessian blocks are to be tested. The default is to test for 1) entire state space, and 2) each dimension individually.
dx	A vector of same length as x. For i'th dimension, nevals values are sampled from a uniform distribution with min/max values equal to x[i]-0.5*dx[i] and x[i]+0.5*dx[i], respectively. Vectors smaller than length(x) are extended as needed by recycling the provided values for dx.
nevals	Number of points in state space, for which validation tests will be performed.
negdef.tol	Lower bound for absolute value of (negative) eigenvalues of Hessian, evaluated at each of the nevals points in the state space. If one or more eigenvalues have absolute values smaller than ngdef.tol, log-density is declared non-log-concave at that point.
...	Other arguments to be passed to fghEval.

Value

sns.check.logdensity returns a list of class sns.check.logdensity, with the following elements:

check.ld.struct	Boolean flag, indicating whether log-density fghEval has one of the 3 forms of output, described above.
numderiv	Integer with values of 0,1,2. A value of 0 means analytical gradient and Hessian have been provided, and thus there is no need for numerical differentiation. 1 means analytical gradient is provided, but Hessian must be calculated numerically. 2 means both gradient and Hessian must be numerically calculated. Users can pass this value to subsequent sns or sns.run calls.
check.length.g	Boolean flag, indicating whether length of gradient vector (element g) returned by fghEval equals length(x).
check.dim.h	Boolean flag, indicating whether number of rows and columns of the Hessian matrix (element h) returned by fghEval equal length(x).
x.mat	Collection of state space vectors (one per row), for which validation tests are performed. It has nevals rows and length(x) columns.
t.evals	Time spent on evaluating fghEval on nevals points chosen randomly in the neighborhood of x, as specified by dx. This includes log-density and, if provided, analytical evaluations of gradient and Hessian.
t.num.evals	Time spent on evaluating the numeric version of fghEval, in which gradient and Hessian are computed numerically, using grad and hessian functions in the numDeriv package. Comparison of this number with t.evals provides the user with insight into the relative speed of numerical differentiation compared to analytical versions.
f.vec	Vector of log-density values for state space vectors listed in x.mat.
g.mat.num	Collection of numerically-computed gradient vectors for state space values listed in x.mat, with the same dimension conventions.
is.g.num.finite	Boolean flag, indicating whether all numerically-computed gradient vectors have finite values.
h.array.num	Collection of numerically-computed Hessian matrices at points listed in x.mat. First dimension is of length nevals, and the remaining two dimensions equal length(x).
is.h.num.finite	Boolean flag, indicating whether all numerically-computed Hessian matrices have finite values.
g.mat	Collection of analytically-computed gradient vectors for state space values listed in x.mat, with the same dimension conventions. This is only available if fghEval has a g field; otherwise NA.
is.g.finite	Boolean flag (if available), indicating whether all analytically-computed gradient vectors have finite values (if available).
g.diff.max	If available, maximum relative difference between analytical and numerical gradient vectors, over all nevals points in x.mat. Relative diference is defined as L2 norm of difference between the two gradient vectors, divided by the L2 norm of the analytical gradient vector.
h.array	If available, collection of analytically-computed Hessian matrices at points listed in x.mat. Dimensional conventions are the same as h.array.num.
is.h.finite	Boolean flag (if available), indicating whether all analytically-computed Hessian matrices have finite values.
h.diff.max	If available, maximum relative difference between analytical and numerical Hessian matrices, over all nevals points in x.mat. Relative difference is defined as the Frobenius norm of difference of analytical and numerical Hessian matrices, divided by the Frobenius norm of analytical Hessian.
is.negdef.num	Boolean flag, indicating whether numerical Hessian is negative-definite at all state space points indicated in x.mat.
is.negdef	Boolean flag, indicating whether analytical Hessian is negative-definite at all state space points indicated in x.mat.

Note

1. Validation tests performed in sns.check.logdensity cannot prove that a log-density is twice-differentiable, or globally concave. However, when e.g. log-density Hessian is seen to be non-negative-definite at one of the points tested, we can definitively say that the Hessian is not globally negative-definite, and therefore sns should not be used for sampling from this distribution. Users must generally consider this function as a supplement to analytical work.

2. See package vignette for more details on SNS theory, software, examples, and performance.

Author(s)

Alireza S. Mahani, Asad Hasan, Marshall Jiang, Mansour T.A. Sharabiani

References

Mahani A.S., Hasan A., Jiang M. & Sharabiani M.T.A. (2016). Stochastic Newton Sampler: The R Package sns. Journal of Statistical Software, Code Snippets, 74(2), 1-33. doi:10.18637/jss.v074.c02

Examples

## Not run: 

# using RegressionFactory for generating log-likelihood and its derivatives
library(RegressionFactory)

loglike.poisson <- function(beta, X, y) {
  regfac.expand.1par(beta, X = X, y = y,
                     fbase1 = fbase1.poisson.log)
}

# simulating data
K <- 5
N <- 1000
X <- matrix(runif(N * K, -0.5, +0.5), ncol = K)
beta <- runif(K, -0.5, +0.5)
y <- rpois(N, exp(X %*% beta))

beta.init <- rep(0.0, K)

my.check <- sns.check.logdensity(beta.init, loglike.poisson
  , X = X, y = y, blocks = list(1:K))
my.check

# mistake in log-likelihood gradient
loglike.poisson.wrong <- function(beta, X, y) {
  ret <- regfac.expand.1par(beta, X = X, y = y,
                            fbase1 = fbase1.poisson.log)
  ret$g <- 1.2 * ret$g
  return (ret)
}
# maximum relative diff in gradient is now much larger
my.check.wrong <- sns.check.logdensity(beta.init
  , loglike.poisson.wrong, X = X, y = y, blocks = list(1:K))
my.check.wrong

# mistake in log-likelihood Hessian
loglike.poisson.wrong.2 <- function(beta, X, y) {
  ret <- regfac.expand.1par(beta, X = X, y = y,
                            fbase1 = fbase1.poisson.log)
  ret$h <- 1.2 * ret$h
  return (ret)
}
# maximum relative diff in Hessian is now much larger
my.check.wrong.2 <- sns.check.logdensity(beta.init
  , loglike.poisson.wrong.2, X = X, y = y, blocks = list(1:K))
my.check.wrong.2


## End(Not run)

sns documentation built on Nov. 2, 2022, 5:15 p.m.