compute_H: Compute all the blocks of H.

In JimSkinner/spca: Structured PCA

Description

Compute all the blocks of H.

Usage

compute_H(X, WHat, muHat, sigSqHat, K)

Arguments

X	Data
WHat	Loadings matrix
K	Prior covariance matrix

Value

H

Examples

set.seed(1)
d=10; k=3; n=1000
X = matrix(rnorm(n*d), ncol=d)
W = matrix(rnorm(d*k), ncol=k)
mu = rnorm(d)
sigSq = rnorm(1)^2
K = cov.SE(matrix(1:10, ncol=1), beta=log(c(2, 3)))

library(numDeriv)
library(Matrix)

#Test that the analytic hessian for mu & sigSq matches numerical Hessian.
H.analytic = stpca:::compute_H(X, W, mu, sigSq, K)
HsigSq.numeric = Matrix(numDeriv::hessian(function(sigSq_) {
  -(log_likelihood(X, W, mu, sigSq_) + log_prior(K, W, sigSq))
}, x=sigSq))
stopifnot(all.equal(H.analytic$sigSq, HsigSq.numeric,
                    tolerance=1e-8))

Hmu.numeric = numDeriv::hessian(function(mu_) {
  -log_likelihood(X, W, mu_, sigSq)
}, x=mu)
stopifnot(isTRUE(all.equal(unname(as.matrix(H.analytic$mu)),
                           Hmu.numeric, tolerance=1e-6)))

JimSkinner/spca documentation built on May 7, 2019, 10:52 a.m.