Nothing
# (C) 2008-2012 Leo Lahti and Olli-Pekka Huovilainen
# All rights reserved.
# FreeBSD License (keep this notice)
# "The mathematics is not there till we put it there."
# - Sir Arthur Eddington, The Philosophy of Physical Science
update.W.singledata <- function (Wt, X, phi, priors = NULL) {
# Optimize W for PFA and PPCA models
# Use cost function to add priors on W
if (is.null(priors)) {
Wvec <- optim(as.vector(Wt), pfa.neg.log.likelihood, X = X, phi = phi, method = "SANN")$par # NelderMead and BFGS do not work as well
} else if (!is.null(priors$W)) {
Wvec <- abs(optim(as.vector(Wt), pfa.cost.regularized, X = X, phi = phi, priors = priors, method = "SANN")$par) # NelderMead and BFGS do not work as well
}
W <- matrix(Wvec, ncol = nrow(X))
}
pfa.cost.regularized <- function (Wvec, phi, X, priors) {
Wvec <- abs(Wvec)
cost.data <- pfa.neg.log.likelihood(Wvec, phi, X)
# -logP for W prior
# wcost <- sum((W$X)^2) * priors$W
# Assuming exponential prior distribution with rate parameter priors$W
#wcost <- 0 # no effect # FIXME: would be faster without this 'if' check
#if ( !is.null(priors$W) && priors$W > 0 ) {
#wcost <- -sum(dexp(Wvec, rate = priors$W, log = TRUE))
#wcost <- wprior.c(Wvec, priors$W)
wcost <- wprior(Wvec, priors$W)
#}
cost.data + wcost
}
wprior <- function (vec, rate) {
-sum(dexp(vec, rate = rate, log = TRUE))
}
#wprior.c <- cmpfun(wprior)
pfa.neg.log.likelihood <- function (Wvec, phi, X) {
# Cost function for W in the PFA model.
# Also applicable for PPCA.
W <- matrix(Wvec, ncol = nrow(X))
# X: features x samples; assuming that this is centered at origo
# X is Y in Rubin-Thayer 1982: this log-likelihood is from Eq. 1 in there
# R <- diag(1, zDimension) # R = I ie. exploratory factor analysis, see Rubin-Thayer Case 1.
# k <- tau2 + t(beta) %*% R %*% beta
# k <- tau2 + t(beta) %*% beta # assuming R = I
# beta <- t(W$total)
# tau2 <- phi$total
# k <- tau2 + t(beta) %*% beta # assuming R = I
# wtw <- W%*%t(W)
# assuming R = I and adding small constant to avoid numerical overflows
k <- t(W)%*%W + phi + diag(1e-18, nrow(phi))
pfacost(ncol(X), k, X)
#pfacost.c(ncol(X), k, X)
}
pfacost <- function (n, k, X) {
-as.numeric(-(n/2)*(determinant(k, logarithm = TRUE)$modulus + sum(diag(cov(t(X)) %*% solve( k )))))
}
cost.W <- function (vec, phi, priors, Dim, Dcov) {
# (C) 2008-2011 Leo Lahti and Olli-Pekka Huovilainen
# All rights reserved.
# FreeBSD License (keep this notice)
# Wx ~ Wy constrained
# no W prior
# Retrieve the actual W and T from the parameter vector
wt <- get.W(vec, Dim)
W <- wt$W
T <- wt$T
# Marginal cost for the whole data set
# integrated over z
# given parameters W, phi
# P(X,Y | W, phi) = integral N(X|Wx*z,phix)*N(Y|Wy*z,phiy)*N(z|0,I)
# We report -logP here
# Data prob. Taken from probCCA paper, section 4, l1
wtw.xy <- W$X%*%t(W$Y)
Sigma <- rbind(cbind(W$X%*%t(W$X) + phi$X*diag(Dim$X), wtw.xy),
cbind(t(wtw.xy),W$Y%*%t(W$Y) + phi$Y*diag(Dim$Y)))
# -logP for the data
cost.data <- log(det(Sigma)) + sum(diag(solve(Sigma)%*%Dcov$total))
# -logP for T prior
tcost <- sum((T - priors$Nm.wxwy.mean)^2) * priors$T.tmp
# -logP for W prior - skip since not used now
# priors$W.tmp <- 1/(2 * Nsamples * priors$W)
# NOTE considerable speed increase in optimize iteration if
# this is calculated outside this function!
#wcost <- sum((W$X)^2) * priors$W.tmp
cost.data + tcost #+ wcost
}
cost.W.exponential <- function (vec, phi, priors = NULL, Dim, Dcov) {
# (C) 2008-2011 Leo Lahti and Olli-Pekka Huovilainen
# All rights reserved.
# FreeBSD License (keep this notice)
# allows exponential prior for W
# in general, Wx != Wy
# remove sign as we assume W always positive here
vec <- abs(vec)
# Retrieve W from the parameter vector
W <- get.W.nonneg(vec, Dim)
# Marginal cost for the whole data set
# integrated over z
# given parameters W, phi
# P(X,Y | W, phi) = integral N(X|Wx*z,phix)*N(Y|Wy*z,phiy)*N(z|0,I)
# report -logP here
# Data prob. Taken from probCCA paper, section 4, l1
wtw.xy <- W$X%*%t(W$Y)
Sigma <- rbind(cbind(W$X%*%t(W$X) + phi$X, wtw.xy),
cbind(t(wtw.xy), W$Y%*%t(W$Y) + phi$Y))
# -logP for the data
det.sigma <- det(Sigma)
if (det.sigma > 1e-6) { # using > 0 caused overflows
cost.data <- log(det.sigma) + sum(diag(solve(Sigma)%*%Dcov$total))
} else {
# In some rare situations det.sigma appears non-positive;
# such solutions are not feasible
cost.data <- 1e300 # 1e309 gives Inf; Inf gives error
}
# -logP for W prior
# wcost <- sum((W$X)^2) * priors$W
# Assuming exponential prior distribution with rate parameter priors$W
wcost <- 0 # no effect # FIXME: would be faster without this 'if' check
if ( !is.null(priors$W) && priors$W > 0 ) {
wcost <- -sum(dexp(vec, rate = priors$W, log = TRUE))
}
cost.data + wcost
}
cost7 <- function (Wvec, phi, Dcov, Dim, priors) {
# (C) 2008-2011 Leo Lahti and Olli-Pekka Huovilainen
# All rights reserved.
# FreeBSD License (keep this notice)
# SimCCA: identical Wx = Wy
# allows W prior
# NOTE: possible to optimize quite much by removing W matrix conversions?
if ( !is.null(priors$W) ) { Wvec <- abs(Wvec) }
W <- get.W.nonneg.identical(Wvec, Dim)
wtw <- W%*%t(W)
Sigma <- rbind(cbind(wtw + phi$X, wtw),
cbind(wtw, wtw + phi$Y))
# Marginal cost for the whole data set
# integrated over z
# given parameters W, phi
# P(X,Y | W, phi) = integral N(X|Wx*z,phix)*N(Y|Wy*z,phiy)*N(z|0,I)
# reporting -logP here
# restrict solutions to cases where det(Sigma)>=0
# -logP for the data
detsigma <- det( Sigma )
if (detsigma > 1e-100) { # crashes with > 0
cost.data <- log(detsigma) + sum(diag(solve(Sigma)%*%Dcov$total))
} else { cost.data <- 1e300 }
# -logP for W prior
# wcost <- sum((W$X)^2) * priors$W
# Assuming exponential prior distribution with rate parameter priors$W
wcost <- 0 # no effect by efault
if ( !is.null(priors$W) ) {
#multiply by 2 to count for both wx and wy
wcost <- -2*sum(dexp(Wvec, rate = priors$W, log = TRUE))
}
cost.data + wcost
}
W.cca.EM <- function (Dcov, M, beta) {
# (C) 2008-2011 Leo Lahti and Olli-Pekka Huovilainen
# All rights reserved.
# FreeBSD License (keep this notice)
# EM update for W, given phi (through M, beta)
# Returns total W, i.e. [Wx; Wy]
#beta <-M%*%t(W$total)%*%phi.inv$total
#M <- solve(t(W)%*%phi.inv%*%W + I) # W$total meant here
ctb <- Dcov$total%*%t(beta)
matrix(ctb%*%solve(M + beta%*%ctb), nrow = nrow(Dcov$total))
}
W.simcca.EM <- function (W, phi, Dim, Dcov) {
# (C) 2008-2011 Leo Lahti and Olli-Pekka Huovilainen
# All rights reserved.
# FreeBSD License (keep this notice)
# CCA update for W, assuming Wx = Wy
# see Bach-Jordan 2005, sec. 4.1 for details
# equations modified from there to match Wx = Wy case
# FIXME: speedup by sharnig M/beta with phi updates as in in W.cca.EM?
# (phi.EM.simcca)
what <- 2*W$X
twp <- t(what)%*%solve(phi$X + phi$Y)
M.w <- solve(twp%*%what + diag(Dim$Z))
beta.w <- M.w%*%twp
ctb <- Dcov$sum%*%t(beta.w)
w <- ctb%*%solve(M.w + beta.w%*%ctb)/2
list(X = w, Y = w, total = rbind(w, w))
}
#######################################
solve.w <- function (Xc, Yc, Cxx, Cyy, dz = NULL) {
# (C) 2008-2011 Leo Lahti and Olli-Pekka Huovilainen
# All rights reserved.
# FreeBSD License (keep this notice)
# FIXME: compare with the other W updates, e.g. W.cca.EM
# assumes Xc, Yc : samples x features, zero-mean features
# Cxx and Cyy are covariances from cov(Xc) and cov(Yc)
# dz shows the desired rank of latent Z
# NOTE: here the dimensions of Xc and Yc do not need to match
# Note: in previous solve.w the input data was features x samples
# Traditional CCA solution (modified from cancor function):
nr <- nrow(Xc)
qx <- qr(Xc)
qy <- qr(Yc)
dx <- qx$rank
dy <- qy$rank
if (dx < ncol(Xc) || dy < ncol(Yc)) {
stop("Unable to calculate pCCA; sample covariance matrix not invertible.")
}
z <- svd(qr.qty(qx, qr.qy(qy, diag(1, nr, dy)))[1L:dx, , drop = FALSE], dx, dy)
xcoef <- backsolve((qx$qr)[1L:dx, 1L:dx, drop = FALSE], z$u)
ycoef <- backsolve((qy$qr)[1L:dy, 1L:dy, drop = FALSE], z$v)
# Solve W using Archambeau06 equations
# Note: only requirement for Q is that Qx%*%t(Qy) = canonical correlations
# Q corresponds to M (dz x dz) in Bach-Jordan 2005, p.8 (before sec 4.1)
Qx <- diag(z$d[1:dz], dz, dz) # dz x dz matrix
#Qy <- diag(1, nrow(Qx)) # also a dz x dz matrix: identity matrix -> omit
# ML estimates for the prob. model W:
dz <- ifelse(is.null(dz), length(z$d), dz)
Wx <- Cxx%*%xcoef[,1:dz]%*%Qx
Wy <- Cyy%*%ycoef[,1:dz]#%*%Qy # Qy is identity matrix -> omit
list(X = Wx, Y = Wy)
}
solve.archambeau <- function (X, Y, Wx, Wy, btb.x, btb.y) {
# (C) 2008-2011 Leo Lahti and Olli-Pekka Huovilainen
# All rights reserved.
# FreeBSD License (keep this notice)
# Use the trick introduced in Archambeau et al., ICML 2006. Robust
# probabilistic projections, and the later correction appendix.
# This is used to retrieve the original principal components from
# the probabilistic CCA model solution i.e. get round the rotational
# invariance problem
# NOTE: Archambeau uses B in a different meaning than Bach-Jordan.
# Here denote his B with B.arch
# btb = B%*%t(B)
Nd <- nrow(Wx)
Bx.arch <- Wx%*%solve(btb.x)%*%t(Wx) + diag(Nd)
By.arch <- Wy%*%solve(btb.y)%*%t(Wy) + diag(Nd)
# R is rotation matrix
R <- eigen((diag(Nd) - solve(Bx.arch))%*%(diag(Nd) - solve(By.arch)))$vector
Ux <- solve(cov(t(X)))%*%Wx%*%solve(matrix.sqrt((diag(Nd)-solve(Bx.arch))))%*%R
Uy <- solve(cov(t(Y)))%*%Wy%*%solve(matrix.sqrt((diag(Nd)-solve(By.arch))))%*%R
list(X = Ux, Y = Uy)
}
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.