Nothing
# X: data matrix ${\bf x}_i, i=1, 2, \ldots, N$
# t: time coordinate
#source("R/spSupport.R")
# if( Sys.info()["sysname"]=="Linux"){
# if(!exists("compileme")) compileme =FALSE
# if(compileme) system("R CMD SHLIB -fopenmp ../src/morph.c")
#
# if( is.loaded("anaForwardC" )) {
# dyn.unload("src/morph.so")
# }
# dyn.load("src/morph.so")
# }else{
# if( is.loaded("anaForwardC" )) {
# dyn.unload("src/gmGeostats.dll")
# }
# dyn.load("src/gmGeostats.dll")
# }
################### auxiliary functions ------------------------------
# anaXt: position of each datum after a time "t": ${\bf x}_i(t)$
anaXt <- function(t,X) {
(1-t)*X
}
# this is a scaling of each row towards the 0-vector:
# anaSigmat: kernel bandwidth at time "t": $\sigma(t)$
anaSigmat <- function(t,sigma0=0.1,sigma1=1) {
(sigma1-sigma0)*t+sigma0
}
# this goes from an initial value to the standard normal distribution
# x: an arbitrary point
# X: one of the "pivots" (data)
# t: time instant
# calculate the normal score of the point w.r.t. the sigmat and Xt
anaScore <- function(x,X,t,...){
(x-anaXt(t,X))/anaSigmat(t,...)
}
exp2 = function(x) exp(-0.5*x^2)
# weights: matrix of weights (rows as many as rows of x, columns as many as rows of Y); by rows sum to one
anaW <- function(x,t,Y,sigma0=sigma0,sigma1=sigma1, funkernel = exp2) {
# sigma = anaSigmat(t,sigma0=sigma0,sigma1=sigma1)
# W<-matrix(exp(-0.5*c((x[rep(1:nrow(x),nrow(Y)),]-anaXt(t,Y)[rep(1:nrow(Y),each=nrow(x)),])/sigma)^2 %*% rep(1,ncol(Y))),nrow=nrow(x))
W<-matrix(funkernel(c(((x[rep(1:nrow(x),nrow(Y)),]-anaXt(t,Y)[rep(1:nrow(Y),each=nrow(x)),])/anaSigmat(t,sigma0=sigma0,sigma1=sigma1))^2%*%rep(1,ncol(Y)))),nrow=nrow(x))
e<- W/c(W%*%rep(1,ncol(W)))
if(any(e)){
warning("sum((e%*%rep(1,ncol(W))-1)^2)<0.0001 not satisfied")
}
e
}
# x: points that are being deformed
# Y: data points
# average speed of deformation of the space at time t
anaV <- function(x,t,Y,sigma0,sigma1=1+sigma0) {
W <- anaW(x,t,Y,sigma0=sigma0,sigma1=sigma1)
i <- rep(1:nrow(x),each=nrow(Y))
j <- rep(1:nrow(Y),nrow(x))
Vind <- matrix(c(t(W))*(-Y[j,]+((sigma1-sigma0)/anaSigmat(t,sigma0=sigma0,sigma1=sigma1))*(x[i,]-anaXt(t,Y)[j,])),nrow=nrow(Y))
Vall <- matrix(rep(1,nrow(Y))%*%Vind,nrow=nrow(x))
Vall
}
################### transformation functions ------------------------------
# # forward
# anaForwardOld <- function(x,Y,sigma0,sigma1=1+sigma0,steps=30,plt=FALSE,sphere=TRUE) {
# if( sphere ) st = sphTrans(Y) else st=function(x,...) x
# x<-st(x)
# Y<-st(Y)
# h=1/steps
# if( plt )
# plot(x[,1:2],pch=".")
# for(i in 1:steps) {
# # Verfahren 2.Ordnung
# v1 = anaV(x,(i-1)*h,Y,sigma0=sigma0,sigma1=sigma1)
# v2 = anaV(x+h*v1,i*h,Y,sigma0=sigma0,sigma1=sigma1)
# x = x+0.5*h*(v1+v2)
# #x = x+ h*anaV(x,i*h,Y,sigma0=sigma0,sigma1=sigma1)
# if(plt)
# points(x[,1:2],pch=".")
# }
# x
# }
#' Forward gaussian anamorphosis
#' forward transformation to multivariate gaussian scores
#'
#' @param x points to be transformed (a matrix)
#' @param Y node points defining the transformation (another matrix, same nr. of columns as `x`)
#' @param sigma0 starting spread of the kernels
#' @param sigma1 final spread of the kernels
#' @param steps number of steps to linearize the transform (default 30 is good)
#' @param plt boolean, do you want to get a plot of the transformation?
#' @param sphere boolean, should the data be pre-Y-spherified first? defaults to true
#' @param weights vector of weights for all computations, length must be equal
#' to number of rows of \code{x}
#'
#' @return a matrix with the gaussian scores; same dimensions of \code{x}
#' @export
#' @seealso [ana()] for defining a function that carries over the transformation
#' (by means of a closure), [anaBackward()] for the explicit back-transformation,
#' [sphTrans()] for defining a function that carries over the spherification of the data
#' @useDynLib gmGeostats anaForwardC
#' @author K. Gerald van den Boogaart, Raimon Tolosana-Delgado
#' @examples
#' data("jura", package="gstat")
#' Y = jura.pred[,c(10,12,13)]
#' plot(compositions::acomp(Y))
#' Ylr = compositions::alr(Y)
#' plot(Ylr)
#' z = anaForward(x=Ylr, Y=Ylr, sigma0=0.1)
#' plot(z, asp=1)
#' shapiro.test(z[,1])
#' shapiro.test(z[,2])
anaForward <- function(x,Y,sigma0,sigma1=1+sigma0,steps=30,plt=FALSE,sphere=TRUE, weights=NULL) {
if(is.null(weights)){
weights = rep(1,nrow(Y))
}
if(length(weights)!=nrow(Y)){
stop("weights provided not compatible with nrow(Y)")
}
# if( sphere ) st = sphTrans(Y) else
st=function(x,...) x
if(is.logical(sphere)){
if( sphere ) st = sphTrans(Y, weights=weights) else st=function(x,...) x
}else if(is.function(sphere)){
st = sphere
}
x<-t(st(x))
Y<-t(st(Y))
#h=1/steps
#if( plt )
# plot(x[,1:2],pch=".")
erg<-.C(anaForwardC,
dimX=checkInt(dim(x),2),
x=checkDouble(x),
dimY=checkInt(dim(Y),2),
y=checkDouble(Y),
wY=checkDouble(weights,ncol(Y)),
steps=checkInt(steps,1),
sigma0=checkDouble(sigma0,1),
sigma1=checkDouble(sigma1,1)#, PACKAGE = "gmGeostats"
)
x<- t(structure(erg$x, dim=dim(x)))
colnames(x) = paste("flow", 1:ncol(x), sep="")
return(x)
}
# anaBackwardOld <- function(x,Y,sigma0,sigma1=1+sigma0,steps=30,plt=FALSE,sphere=TRUE) {
# if( sphere ) st = sphTrans(Y) else st=function(x,...) x
# Y<-st(Y)
# h=1/steps
# if( plt )
# plot(x[,1:2],pch=".")
# for(i in 1:steps) {
# # Verfahren 2.Ordnung
# v1 = anaV(x,1-(i-1)*h,Y,sigma0=sigma0,sigma1=sigma1)
# v2 = anaV(x-h*v1,1-i*h,Y,sigma0=sigma0,sigma1=sigma1)
# x = x-0.5*h*(v1+v2)
# # x = x- h*anaV(x,1-i*h,Y,sigma0=sigma0,sigma1=sigma1)
# if(plt)
# points(x[,1:2],pch=".")
# }
# st(x,inv=TRUE)
# }
#' Backward gaussian anamorphosis
#' backward transformation to multivariate gaussian scores
#'
#' @param x matrix of gaussian scores to be back-transformed
#' @param Y node points defining the transformation (a matrix, same nr of columns)
#' @param sigma0 starting spread of the kernels in the forward transform
#' @param sigma1 final spread of the kernels in the forward transform
#' @param steps number of steps to linearize the transform (default 30 is good)
#' @param plt boolean, do you want to get a plot of the transformation?
#' @param sphere boolean, should the data be taken as pre-Y-spherified? defaults to true
#' @param weights vector of weights for all computations, length must be equal
#' to number of rows of \code{x}
#'
#' @return a matrix with the scores back-transformed to the same scale as \code{Y}; same dimensions of \code{x}
#' @seealso [ana()] for defining a function that carries over the transformation
#' (by means of a closure), [anaBackward()] for the explicit back-transformation,
#' [sphTrans()] for defining a function that carries over the spherification of the data
#' @export
#' @useDynLib gmGeostats anaBackwardC
#' @author K. Gerald van den Boogaart, Raimon Tolosana-Delgado
#' @examples
#' data("jura", package="gstat")
#' Y = jura.pred[,c(10,12,13)]
#' plot(compositions::acomp(Y))
#' Ylr = compositions::alr(Y)
#' Xns = matrix(rnorm(500), ncol=2)
#' plot(Ylr)
#' points(Xns, col=2, pch=4)
#' Xlr = anaBackward(x=Xns, Y=Ylr, sigma0=0.1)
#' qqplot(Xlr[,1], Ylr[,1])
#' qqplot(Xlr[,2], Ylr[,2])
#' qqplot(Xlr[,1]+Xlr[,2], Ylr[,1]+Ylr[,2])
anaBackward <- function(x,Y,sigma0,sigma1=1+sigma0,steps=30,plt=FALSE,sphere=TRUE, weights=NULL) {
if(is.null(weights)){
weights = rep(1,nrow(Y))
}
if(length(weights)!=nrow(Y)){
stop("weights provided not compatible with nrow(Y)")
}
#if( sphere ) st = sphTrans(Y) else
st=function(x,...) x
if(is.logical(sphere)){
if( sphere ) st = sphTrans(Y, weights=weights) else st=function(x,...) x
}else if(is.function(sphere)){
st = sphere
}
Y<-t(st(Y))
x <- t(x)
erg<-.C(anaBackwardC,
dimX=checkInt(dim(x),2),
x=checkDouble(x),
dimY=checkInt(dim(Y),2),
y=checkDouble(Y),
wY=checkDouble(weights,ncol(Y)),
steps=checkInt(steps,1),
sigma0=checkDouble(sigma0,1),
sigma1=checkDouble(sigma1,1) #, PACKAGE = "gmGeostats"
)
x<- t(structure(erg$x, dim=dim(x)))
st(x,inv=TRUE)
}
# pre-spherification of the data (through SVD)
#' Spherifying transform
#' Compute a transformation that spherifies a certain data set
#'
#' @param Y data set defining the spherifization
#' @param ... extra arguments for generic functionality
#'
#' @return a function with arguments \code{(x, inv=FALSE)}, where \code{x} will be the
#' data to apply the transformation to, and \code{inv=FALSE} will indicate if the direct
#' or the inverse transformation is desired.
#' This function applied to the same data returns a translated, rotated and scaled, so that
#' the new scores are centered, have variance 1, and no correlation.
#' @export
#' @importFrom utils methods
#' @importFrom stats cov.wt
#'
#' @aliases sphTrans.default
#' @seealso ana, anaBackward, sphTrans
#' @author K. Gerald van den Boogaart, Raimon Tolosana-Delgado
#'
#' @examples
#' library(compositions)
#' data("jura", package="gstat")
#' Y = acomp(jura.pred[,c(10,12,13)])
#' oldpar = par(mfrow = c(1,1))
#' plot(Y)
#' sph = sphTrans(Y)
#' class(sph)
#' z = sph(Y)
#' plot(z)
#' par(oldpar)
#' cor(cbind(z, ilr(Y)))
#' colMeans(cbind(z, ilr(Y)))
sphTrans <- function(Y,...) UseMethod("sphTrans",Y)
#' @describeIn sphTrans Spherifying transform
#' @method sphTrans default
#' @export
#' @param weights weights to incorporate in the compuations, length=nrow(Y)
#' @param p dimensions to be considered structural (useful for filtering noise)
sphTrans.default <- function(Y, weights=NULL, p=1:ncol(Y), ...) {
cY = class(Y)[1]
existsIdt = sub("idt.", "", as.character(methods("idt")))
if(is.null(weights)){
if("data.frame" %in% cY){
b <- gmApply(Y,2,mean)
}else if(length(intersect(cY, existsIdt))>0){
b <- mean(idt(Y))
}else{
b <- gmApply(Y,2,mean)
}
vY = var(idt(Y))
}else if(length(weights)==nrow(Y)){
aux = cov.wt(idt(Y),wt=weights)
b = aux$center
vY = aux$cov
}else{
stop("meaningless weights provided, wrong length?")
}
p = p[(p<=ncol(vY))]
SVD <- svd(vY)
A <- with(SVD, v[,p] %*% diag(1/sqrt(d[p])) %*% t(u[,p]))
Ai <- with(SVD, u[,p] %*% diag(sqrt(d[p])) %*% t(v[,p]))
f<-function(x,inv=FALSE) {
if(isTRUE(inv)){
if(class(x)[1]=="matrix"){
out = t((Ai%*%t(x))+unclass(b))
}else{
out = Ai%*%x+b
}
return(idtInv(out, orig=Y))
}else{
if(class(x)[1]=="matrix"){
out = t(A%*%(t(x)-unclass(b)))
}else{
out = unclass(A %*% (idt(x)-b))
attr(out, "orig") = NULL
attr(out, "V") = NULL
}
return(out)
}
}
return(f)
}
#' Flow anamorphosis transform
#' Compute a transformation that gaussianizes a certain data set
#'
#' @param Y data set defining the gaussianization
#' @param sigma0 starting spread of the kernels
#' @param sigma1 final spread of the kernels
#' @param steps number of steps to linearize the transform (default 30 is good)
#' @param sphere boolean, should the transform include a spherifization step,
#' making \code{Y} spherical?
#' @param weights weights to incorporate in the compuations, length=nrow(Y)
#'
#' @return a function with arguments \code{(x, inv=FALSE)}, where \code{x} will be the
#' data to apply the transformation to, and \code{inv=FALSE} will indicate if the direct
#' or the inverse transformation is desired
#' @export
#' @seealso anaForward, anaBackward, sphTrans
#' @author K. Gerald van den Boogaart
#' @examples
#' library(compositions)
#' data("jura", package="gstat")
#' Y = acomp(jura.pred[,c(10,12,13)])
#' plot(Y)
#' anafun = ana(Y)
#' class(anafun)
#' z = anafun(Y)
#' plot(z)
#' y = anafun(z, inv=TRUE)
#' plot(data.frame(orig=Y,recalc=y))
ana <-function(Y,sigma0=0.1,sigma1=1,steps=30,sphere=TRUE, weights=NULL) {
f<-function(x,inv=FALSE,...){
if( isTRUE(inv) ) {
anaBackward(x,Y,sigma0=sigma0,sigma1=sigma1,sphere=sphere,weights=weights,...)
} else {
anaForward(x,Y,sigma0=sigma0,sigma1=sigma1,sphere=sphere,weights=weights,...)
}
}
f
}
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.