jranke / mkin
1
# This file is part of the R package mkin
2

3
# mkin is free software: you can redistribute it and/or modify it under the
4
# terms of the GNU General Public License as published by the Free Software
5
# Foundation, either version 3 of the License, or (at your option) any later
6
# version.
7

8
# This program is distributed in the hope that it will be useful, but WITHOUT
9
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
10
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
11
# details.
12

13
# You should have received a copy of the GNU General Public License along with
14
# this program. If not, see <http://www.gnu.org/licenses/>
15

16
#' Function to perform isometric log-ratio transformation
17
#' 
18
#' This implementation is a special case of the class of isometric log-ratio
19
#' transformations.
20
#' 
21
#' @aliases ilr invilr
22
#' @param x A numeric vector. Naturally, the forward transformation is only
23
#'   sensible for vectors with all elements being greater than zero.
24
#' @return The result of the forward or backward transformation. The returned
25
#'   components always sum to 1 for the case of the inverse log-ratio
26
#'   transformation.
27
#' @author René Lehmann and Johannes Ranke
28
#' @seealso Another implementation can be found in R package
29
#'   \code{robCompositions}.
30
#' @references Peter Filzmoser, Karel Hron (2008) Outlier Detection for
31
#'   Compositional Data Using Robust Methods. Math Geosci 40 233-248
32
#' @keywords manip
33
#' @examples
34
#' 
35
#' # Order matters
36
#' ilr(c(0.1, 1, 10))
37
#' ilr(c(10, 1, 0.1))
38
#' # Equal entries give ilr transformations with zeros as elements
39
#' ilr(c(3, 3, 3))
40
#' # Almost equal entries give small numbers
41
#' ilr(c(0.3, 0.4, 0.3))
42
#' # Only the ratio between the numbers counts, not their sum
43
#' invilr(ilr(c(0.7, 0.29, 0.01)))
44
#' invilr(ilr(2.1 * c(0.7, 0.29, 0.01)))
45
#' # Inverse transformation of larger numbers gives unequal elements
46
#' invilr(-10)
47
#' invilr(c(-10, 0))
48
#' # The sum of the elements of the inverse ilr is 1
49
#' sum(invilr(c(-10, 0)))
50
#' # This is why we do not need all elements of the inverse transformation to go back:
51
#' a <- c(0.1, 0.3, 0.5)
52
#' b <- invilr(a)
53
#' length(b) # Four elements
54
#' ilr(c(b[1:3], 1 - sum(b[1:3]))) # Gives c(0.1, 0.3, 0.5)
55
#' 
56
#' @export
57
ilr <- function(x) {
58 1
  z <- vector()
59 1
  for (i in 1:(length(x) - 1)) {
60 1
    z[i] <- sqrt(i/(i+1)) * log((prod(x[1:i]))^(1/i) / x[i+1])
61
  }
62 1
  return(z)
63
}
64

65
#' @rdname ilr
66
#' @export
67
invilr<-function(x) {
68 1
  D <- length(x) + 1
69 1
  z <- c(x, 0)
70 1
  y <- rep(0, D)
71 1
  s <- sqrt(1:D*2:(D+1))
72 1
  q <- z/s
73 1
  y[1] <- sum(q[1:D])
74 1
  for (i in 2:D) {
75 1
    y[i] <- sum(q[i:D]) - sqrt((i-1)/i) * z[i-1]
76
  }
77 1
  z <- vector()
78 1
  for (i in 1:D) {
79 1
    z[i] <- exp(y[i])/sum(exp(y))
80
  }
81

82
  # Work around a numerical problem with NaN values returned by the above
83
  # Only works if there is only one NaN value: replace it with 1
84
  # if the sum of the other components is < 1e-10
85 1
  if (sum(is.na(z)) == 1 && sum(z[!is.na(z)]) < 1e-10)
86 0
    z = ifelse(is.na(z), 1, z)
87

88 1
  return(z)
89
}

Read our documentation on viewing source code .

Loading