R/RDD-OOK_AuxiliaryFunctions-PDFs.R

In alexdidkovskyi/RDDOOK3: Random Domain Decompositions for Object-Oriented Kriging

############################################################
# Random Domain Decompositions for object-oriented Kriging #
#                 over complex domains                     #
#       (A. Menafoglio, G. Gaetani, P. Secchi)             #
#      Supplementary material submitted to JCGS            #
#               Auxiliary functions                        #
############################################################

# This code contains the auxiliary functions needed to work with
# PDF-data based on the Bayes space geometry. 
# Authors: A. Menafoglio, P. Secchi, Politecnico di Milano
# Contact: alessandra.menafoglio@polimi.it
# Last update: 29th Nov. 2017

####

trapzc <- function(t_step, y){ 
  # compute the integral of y with step "t_step"

  return(t_step*(y[1]/2 + sum(y[2:(length(y)-1)]) + y[length(y)]/2 ))
}

trapzc(c(1,2),0)

clr <- function(density, z, z_step){ 
  # transform a density to a clr

  return(log(density) - trapzc(z_step, log(density))/(max(z)-min(z)))
}

clr.mvt <- function(density.df, z, z_step){ 
  # transform a dataset of densities to clr

  if(length(z) != dim(density.df)[1]) {
    N_samples <- dim(density.df)[1]
    n <- dim(density.df)[2]
    density.df <- t(density.df)
  } else{
    N_samples <- dim(density.df)[2]
    n <- dim(density.df)[1]  
  }
  
  
  #res <- matrix( unlist(lapply(1:N_samples,function(x) log(density.df[,x])-trapzc(z_step,log(density.df[,x]))/(max(z)-min(z)))), ncol = N_samples, nrow =n)
  res <- matrix(NA, ncol = N_samples, nrow = n)
  for(i in 1:N_samples)
    res[,i] <- log(density.df[,i])-trapzc(z_step,log(density.df[,i]))/(max(z)-min(z))
  return(res)
}

clr2density <- function(clr, z, z_step){
  # back-transform a clr to a density

  if(is.fd(clr))
    return(exp(eval.fd(z, clr))/trapzc(z_step, exp(eval.fd(z, clr))))
  else
    return(exp(clr)/trapzc(z_step,exp(clr)))
}

clr2density.mv <- function(clr.df, z,z_step){
  # transform a dataset of clr into densities
  
  N_samples <- dim(clr.df)[2]
  dens.df <- matrix(0, nrow = length(z), ncol = N_samples)
  
  for(j in 1:N_samples) {
    dens.df[,j] <- clr2density(clr.df[,j], z, z_step)
  }
  return(dens.df)
}


mean.nr <- function(density.df, z, z_step){ 
  # Compute the mean of the distribution given its density
  
  N_samples <- dim(density.df)[2]
  mean.df <- rep(0, N_samples)
  
  for(j in 1:N_samples){
    mean.df[j] <- trapzc(z_step, z*density.df[,j])
  }
  return(mean.df)
}

var.nr <- function(density.df, z, z_step ){
  # Compute the variance of the distribution given its density
  
  N_samples <- dim(density.df)[2]
  var.df <- rep(0, N_samples)
  mean.df <- rep(0, N_samples)
  
  for(j in 1:N_samples){
    mean.df[j] <- trapzc(z_step, z*density.df[,j])
    var.df[j] <- trapzc(z_step, z^2*density.df[,j]) - mean.df[j]^2
  }
  return(var.df)
}

quantile.dens <- function(t, density, p){ 
  # Compute the quantile of a given order of the distribution from its density
  
  t_step <- t[2] - t[1]
  cdf <- rep(0, length(t)-1)
  for(i in 2:length(t))
    cdf[i] <- trapzc(t_step, density[1:i])
  
  id.min <- max(which(cdf < p))
  t.min <- t[id.min]
  id.max <- min(which(cdf > p))
  t.max <- t[id.max]
  
  w <- (cdf[id.max] - p) / (cdf[id.max] - cdf[id.min])
  q <- t.min*(1-w) + t.max*w
  
  return(q)
}

quantile.mv <- function(t, density, p){ 
  # Compute the quantile of given order of a set 
  # of distributions from their densities
  
  if(length(t) != dim(density)[1])
    density <- t(density)
  q <- rep(NA, dim(density)[2])
  for(j in 1:dim(density)[2]) {
    q[j] <- quantile.dens(t, density[,j], p)
  }
  return(q)
}

p.dens <- function(t, density, t.f){ 
  # Compute the probability on the left of t.f given the density
  
  t_step <- t[2] - t[1]
  id.max <- max(which(t < t.f))
  pf <- trapzc(t_step, density[1:id.max])
  pf1 <- trapzc(t_step, density[1:(id.max + 1)])
  deltapf <- pf1 - pf
  pf <- pf + deltapf*(t.f - t[id.max])/(t_step)
}

p.dens.mv <- function(t, density, t.f){ 
  # Compute the probabilities on the left of t.f from a set of densities
  
  if(length(t) != dim(density)[1])
    density <- t(density)
  p <- rep(NA, dim(density)[2])
  for(j in 1:dim(density)[2]){
    p[j] <-  p.dens(t, density[,j], t.f)
  }
  return(p)
}