initial commit

DESCRIPTION

+Package: ivegGaussr
+Title: Distribution of Gaussian ratios, for the study of vegetation indices
+Version: 0.0.0.9000
+Authors@R: c(person("Pierre", "Santagostini", role = c("aut", "cre"), email = "pierre.santagostini@institut-agro.fr"),
+    person("Angélina", "El Ghaziri", role = "aut", email = "angelina.elghaziri@institut-agro.fr"),
+    person("Nizar", "Bouhlel", "nizar.bouhlel@institut-agro.fr", role = "aut")
+    )
+Maintainer: Pierre Santagostini <pierre.santagostini@institut-agro.fr>
+Description: Tools for the manipulation of Gaussian ratios: parameter estimation, simulation. These can be applied to the study of vegetation indices.
+License: GPL (>= 3)
+Encoding: UTF-8
+Roxygen: list(markdown = TRUE)
+RoxygenNote: 7.3.2

NAMESPACE

+# Generated by roxygen2: do not edit by hand
+
+export(EM)
+export(dz)
+export(kummerM)
+export(lnpochhammer)
+export(pochhammer)
+importFrom(graphics,plot)

R/EM.R

+EM <- function(z, eps = 1e-06) {
+  #' Estimation of the Distribution Parameters of the Ratio of two Gaussian Distributions
+  #'
+  #' Estimation of the parameters of the distribution of a ratio of two distributions \eqn{Z = \frac{X}{Y}},
+  #' \eqn{X} and \eqn{Y} being distributed according to Gaussian distributions,
+  #' using the EM (estimation-maximization) algorithm.
+  #'
+  #' @aliases EM
+  #'
+  #' @usage EM(z, eps = 1e-6)
+  #' @param z numeric matrix or data frame.
+  #' @param eps numeric. Precision for the estimation of the parameters.
+  #' @return A list of 5 elements:
+  #' \itemize{
+  #' \item \code{x}: the mean \eqn{\hat{\mu}_x} and standard deviation \eqn{\hat{\sigma}_x} of the first distribution.
+  #' \item \code{y}: the mean \eqn{\hat{\mu}_y} and standard deviation \eqn{\hat{\sigma}_y} of the first distribution.
+  #' \item \code{beta}, \code{rho}, \code{delta}: the parameters of the \eqn{Z} distribution:
+  #' \eqn{\hat{\delta}_y = \frac{\hat{\sigma}_y}{\hat{\mu}_y}},
+  #' \eqn{\hat{\beta} = \frac{\hat{\mu}_x}{\hat{\mu}_y}},
+  #' \eqn{\hat{\rho} = \frac{\hat{\sigma}_y}{\hat{\sigma}_x}}.
+  #' }
+  #' with two attributes \code{attr(, "epsilon")} (precision of the result) and \code{attr(, "k")} (number of iterations).
+  #'
+  #' @details The  parameters \eqn{\beta} are computed
+  #' using the method presented in Pascal et al., using an iterative algorithm.
+  #'
+  #' The precision for the estimation of \code{beta} is given by the \code{eps} parameter.
+  #'
+  #' @author Pierre Santagostini, Nizar Bouhlel
+  #'
+  #' @references  F. Pascal, L. Bombrun, J.Y. Tourneret, Y. Berthoumieu. Parameter Estimation For Multivariate Generalized Gaussian Distribution.
+  #' IEEE Trans. Signal Processing, vol. 61 no. 23, p. 5960-5971, Dec. 2013.
+  #' \doi{DOI:10.1109/TSP.2013.2282909}
+  #'
+  #' @seealso \code{\link{dz}}: probability density of a Gaussian ratio.
+
+  #' @importFrom graphics plot
+  #' @export
+
+  kummA <- function(x) {
+    Re(kummerM(x, 2, 1.5))/Re(kummerM(x, 1, 0.5))
+  }
+
+  kummB <- function(x) {
+    Re(kummerM(x, 2, 0.5))/Re(kummerM(x, 1, 0.5))
+  }
+
+  # Number of observations
+  n <- length(z)
+
+  #initialisation
+  mux <- 60; muy <- 400
+  variancex <- 20^2; variancey <- 85^2
+  iter <- 0
+  diff <- 100
+
+  while (diff > eps) {
+    iter <- iter + 1
+
+    mu <- 0.5*((z^2/variancex[iter])+1/variancey[iter])
+    gam <- muy[iter]/variancey[iter]+(mux[iter]/variancex[iter])*z
+
+    A <- sapply(gam^2/(4*mu), kummA)
+
+    B <- sapply(gam^2/(4*mu), kummB)
+
+    mux[iter+1] <- (1/n)*sum(z*(gam/mu)*A)
+    variancex[iter+1] <- (1/n)*sum((z^2)*(1/mu)*B)-mux[iter+1]^2
+
+    muy[iter+1] <- (1/n)*sum((gam/mu)*A)
+    variancey[iter+1] <- (1/n)*sum((1/mu)*B )-muy[iter+1]^2
+
+    diff <- max(abs(c(mux[iter+1]-mux[iter],variancex[iter+1]-variancex[iter],muy[iter+1]-muy[iter],variancey[iter+1]-variancey[iter])))
+  }
+
+  sigmax <- sqrt(variancex); sigmay <- sqrt(variancey)
+
+  res <- list()
+  res$x <- list(mu = mux, sigma = sigmax)
+  res$y <- list(mu = muy, sigma = sigmay)
+
+  res$beta <- mux/muy
+  res$rho <- sigmay/sigmax
+  res$delta <- sigmay/muy
+
+  attr(res, "k") <- iter
+  attr(res, "epsilon") <- eps
+
+  return(res)
+}

R/dz.R

+dz <- function (z, bet, rho, delta) {
+  #' Ratio of two Gaussian distributions
+  #'
+  #' Density of the ratio of two Gaussian distributions.
+  #'
+  #' @aliases dz
+  #'
+  #' @usage dz(z, bet, rho, delta)
+  #' @param z length \eqn{p} numeric vector.
+  #' @param bet,rho,delta numeric values. The parameters \eqn{(\beta, \rho, \delta)} of the distribution, see Details.
+  #'
+  #' @details Let two random variables
+  #' \eqn{X \sim N(\mu_x, \sigma_x)} and \eqn{Y \sim N(\mu_y, \sigma_y)}.
+  #' If we denote \eqn{\beta = \frac{\mu_x}{\mu_y}}, \eqn{\rho = \frac{\sigma_y}{\sigma_x}} and \eqn{\delta_y = \frac{\sigma_y}{\mu_y}},
+  #' the probability distribution function of the ratio \eqn{Z = \frac{X}{Y}}
+  #' is given by: \deqn{\displaystyle{ f_Z(z; \beta, \rho, \delta_y) = \frac{\rho}{\pi (1 + \rho^2 z^2)} \left[ \exp{\left(-\frac{\rho^2 \beta^2 + 1}{2\delta_y^2}\right)} + \sqrt{\frac{\pi}{2}} \ q \ \text{erf}\left(\frac{q}{\sqrt{2}}\right) \exp\left(-\frac{\rho^2 (z-\beta)^2}{2 \delta_y^2 (1 + \rho^2 z^2)}\right) \right] }}
+  #' with \eqn{ q = \frac{1 + \beta \rho^2 z}{\delta_y \sqrt{1 + \rho^2 z^2}} }
+  #' and \eqn{ \text{erf}\left(\frac{q}{\sqrt{2}}\right) = \frac{2}{\sqrt{\pi}} \int_0^{\frac{q}{\sqrt{2}}}{\exp{(-t^2)}\ dt} }
+  #'
+  #' @author Pierre Santagostini, Nizar Bouhlel
+  #'
+  #' @references El Ghaziri, A., Bouhlel, N., Sapoukhina, N., Rousseau, D.,
+  #' On the importance of non-Gaussianity in chlorophyll fluorescence imaging.
+  #' Remote Sensing 15(2), 528 (2023).
+  #' \doi{10.3390/rs15020528}
+  #'
+  #' Díaz-Francés, E., Rubio, F.J.,
+  #' On the existence of a normal approximation to the distribution of the ratio of two independent normal random variables.
+  #' Stat Papers 54, 309–323 (2013).
+  #' \doi{10.1080/03610929808832115}
+  #'
+  #' @examples
+  #' # First example
+  #' beta1 <- 0.15
+  #' rho1 <- 5.75
+  #' delta1 <- 0.22
+  #' dz(0, bet = beta1, rho = rho1, delta = delta1)
+  #' dz(0.5, bet = beta1, rho = rho1, delta = delta1)
+  #' curve(dz(x, bet = beta1, rho = rho1, delta = delta1), from = -0.1, to = 0.7)
+  #'
+  #' # Second example
+  #' beta2 <- 0.24
+  #' rho2 <- 4.21
+  #' delta2 <- 0.25
+  #' dz(0, bet = beta2, rho = rho2, delta = delta2)
+  #' dz(0.5, bet = beta2, rho = rho2, delta = delta2)
+  #' curve(dz(x, bet = beta2, rho = rho2, delta = delta2), from = -0.1, to = 0.7)
+  #'
+  #' @export
+
+  q <- (1+bet*rho^2*z)/(delta*sqrt(1+rho^2*z^2))
+  denom = rho/(pi*(1+rho^2*z^2))*(
+    exp(-(rho^2*bet^2+1)/(2*delta^2)) +
+      sqrt(pi/2)*q*.erf(q/sqrt(2))*exp(-0.5*(rho^2*(z-bet)^2)/(delta^2*(1+rho^2*z^2)))
+  )
+  return(denom)
+}
+
+.erf <- Vectorize(function(x) 2*pnorm(sqrt(2)*x) - 1)

R/kummerM.R

+kummerM <- function(a, b, z, eps = 1e-06) {
+  #' Lauricella \eqn{D}-Hypergeometric Function
+  #'
+  #' Computes the Kummer's function, or confluent hypergeometric function.
+  #'
+  #' @aliases lauricella
+  #'
+  #' @usage kummerM(a, b, z, eps = 1e-06)
+  #' @param a numeric.
+  #' @param b numeric
+  #' @param z numeric.
+  #' @param eps numeric. Precision for the sum (default 1e-06).
+  #' @return A numeric value: the value of the Kummer's function,
+  #' with two attributes \code{attr(, "epsilon")} (precision of the result) and \code{attr(, "k")} (number of iterations).
+  #'
+  #' @details The Kummer's confluent hypergeometric function is given by:
+  #' \deqn{\displaystyle{_1 F_1\left(a, b; z\right) = \sum_{n = 0}^{+\infty}{ \frac{ (a)_n }{ (b)_n } \frac{z^n}{n!} }}}
+  #'
+  #' where \eqn{(z)_p} is the Pochhammer symbol (see \code{\link{pochhammer}}).
+  #'
+  #' The \code{eps} argument gives the required precision for its computation.
+  #' It is the \code{attr(, "epsilon")} attribute of the returned value.
+  #'
+  #' @author Pierre Santagostini, Nizar Bouhlel
+  #' 
+  #' N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions.
+  #' IEEE Signal Processing Letters, vol. 30, pp. 1672-1676, October 2023.
+  #' \doi{10.1109/LSP.2023.3324594}
+  #' @export
+  
+  d <- 1
+  
+  n <- 0
+  res <- 0
+  
+  while ((max(d) > eps) & (!is.nan(d))) {
+    res <- res + d
+    n <- n + 1
+    d <- exp(sapply(z, function(x)
+      Re(lnpochhammer(a, n)*n*log(x) - (lnpochhammer(b, n) + lfactorial(n)))
+    ))
+  }
+  
+  attr(res, "k") <- n
+  attr(res, "epsilon") <- d
+  return(res)
+}
\ No newline at end of file

R/lnpochhammer.R

+lnpochhammer <- function(x, n) {
+  #' Logarithm of the Pochhammer Symbol
+  #'
+  #' Computes the logarithm of the Pochhammer symbol.
+  #'
+  #' @aliases lnpochhammer
+  #'
+  #' @usage lnpochhammer(x, n)
+  #' @param x numeric.
+  #' @param n positive integer.
+  #' @return Numeric value. The logarithm of the Pochhammer symbol.
+  #' @details The Pochhammer symbol is given by:
+  #' \deqn{ \displaystyle{ (x)_n = \frac{\Gamma(x+n)}{\Gamma(x)} = x (x+1) ... (x+n-1) } }
+  #' So, if \eqn{n > 0}:
+  #' \deqn{ \displaystyle{ log\left((x)_n\right) = log(x) + log(x+1) + ... + log(x+n-1) } }
+  #' 
+  #' If \eqn{n = 0}, \eqn{\displaystyle{ log\left((x)_n\right) = log(1) = 0}}
+  #' @seealso \code{\link{pochhammer}},
+  #' \code{\link{lauricella}}
+  #' @author Pierre Santagostini, Nizar Bouhlel
+  #'
+  #' @examples
+  #' lnpochhammer(2, 0)
+  #' lnpochhammer(2, 1)
+  #' lnpochhammer(2, 3)
+  #'
+  #' @export
+  
+  # Arguments:
+  #   - x: numeric
+  #   - n: integer
+  
+  if (n < 0) {
+    stop("n must be non negative")
+  }
+  
+  if (n == 0) {
+    return(0)
+  }
+  
+  if (n > 0) {
+    y <- x + (1:n) - 1
+    return(sum(log(abs(y)) + (y < 0)*pi*1i))
+    # y <- matrix(x, nrow = n, ncol = length(x), byrow = TRUE) +
+    #   matrix(1:n, nrow = n, ncol = length(x), byrow = FALSE) - 1
+    # return(colSums(log(abs(y)) + (y < 0)*pi*1i))
+  }
+}

R/pochhammer.R

+pochhammer <- function(x, n) {
+  #' Pochhammer Symbol
+  #'
+  #' Computes the Pochhammer symbol.
+  #'
+  #' @aliases pochhammer
+  #'
+  #' @usage pochhammer(x, n)
+  #' @param x numeric.
+  #' @param n positive integer.
+  #' @return Numeric value. The value of the Pochhammer symbol.
+  #' @details The Pochhammer symbol is given by:
+  #' \deqn{ \displaystyle{ (x)_n = \frac{\Gamma(x+n)}{\Gamma(x)} = x (x+1) ... (x+n-1) } }
+  #' @seealso \code{\link{lauricella}}
+  #' @author Pierre Santagostini, Nizar Bouhlel
+  #'
+  #' @examples
+  #' pochhammer(2, 0)
+  #' pochhammer(2, 1)
+  #' pochhammer(2, 3)
+  #'
+  #' @export
+
+  # Arguments:
+  #   - x: numeric
+  #   - n: integer
+
+  if (n < 0) {
+    stop("n must be non negative")
+  }
+
+  if (n == 0) {
+    return(1)
+  }
+
+  if (n > 0) {
+    return(prod(x + (1:n) - 1))
+  }
+}

man/EM.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/EM.R
+\name{EM}
+\alias{EM}
+\title{Estimation of the Distribution Parameters of the Ratio of two Gaussian Distributions}
+\usage{
+EM(z, eps = 1e-6)
+}
+\arguments{
+\item{z}{numeric matrix or data frame.}
+
+\item{eps}{numeric. Precision for the estimation of the parameters.}
+}
+\value{
+A list of 5 elements:
+\itemize{
+\item \code{x}: the mean \eqn{\hat{\mu}_x} and standard deviation \eqn{\hat{\sigma}_x} of the first distribution.
+\item \code{y}: the mean \eqn{\hat{\mu}_y} and standard deviation \eqn{\hat{\sigma}_y} of the first distribution.
+\item \code{beta}, \code{rho}, \code{delta}: the parameters of the \eqn{Z} distribution:
+\eqn{\hat{\delta}_y = \frac{\hat{\sigma}_y}{\hat{\mu}_y}},
+\eqn{\hat{\beta} = \frac{\hat{\mu}_x}{\hat{\mu}_y}},
+\eqn{\hat{\rho} = \frac{\hat{\sigma}_y}{\hat{\sigma}_x}}.
+}
+with two attributes \code{attr(, "epsilon")} (precision of the result) and \code{attr(, "k")} (number of iterations).
+}
+\description{
+Estimation of the parameters of the distribution of a ratio of two distributions \eqn{Z = \frac{X}{Y}},
+\eqn{X} and \eqn{Y} being distributed according to Gaussian distributions,
+using the EM (estimation-maximization) algorithm.
+}
+\details{
+The  parameters \eqn{\beta} are computed
+using the method presented in Pascal et al., using an iterative algorithm.
+
+The precision for the estimation of \code{beta} is given by the \code{eps} parameter.
+}
+\references{
+F. Pascal, L. Bombrun, J.Y. Tourneret, Y. Berthoumieu. Parameter Estimation For Multivariate Generalized Gaussian Distribution.
+IEEE Trans. Signal Processing, vol. 61 no. 23, p. 5960-5971, Dec. 2013.
+\doi{DOI:10.1109/TSP.2013.2282909}
+}
+\seealso{
+\code{\link{dz}}: probability density of a Gaussian ratio.
+}
+\author{
+Pierre Santagostini, Nizar Bouhlel
+}

man/dz.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/dz.R
+\name{dz}
+\alias{dz}
+\title{Ratio of two Gaussian distributions}
+\usage{
+dz(z, bet, rho, delta)
+}
+\arguments{
+\item{z}{length \eqn{p} numeric vector.}
+
+\item{bet, rho, delta}{numeric values. The parameters \eqn{(\beta, \rho, \delta)} of the distribution, see Details.}
+}
+\description{
+Density of the ratio of two Gaussian distributions.
+}
+\details{
+Let two random variables
+\eqn{X \sim N(\mu_x, \sigma_x)} and \eqn{Y \sim N(\mu_y, \sigma_y)}.
+If we denote \eqn{\beta = \frac{\mu_x}{\mu_y}}, \eqn{\rho = \frac{\sigma_y}{\sigma_x}} and \eqn{\delta_y = \frac{\sigma_y}{\mu_y}},
+the probability distribution function of the ratio \eqn{Z = \frac{X}{Y}}
+is given by: \deqn{\displaystyle{ f_Z(z; \beta, \rho, \delta_y) = \frac{\rho}{\pi (1 + \rho^2 z^2)} \left[ \exp{\left(-\frac{\rho^2 \beta^2 + 1}{2\delta_y^2}\right)} + \sqrt{\frac{\pi}{2}} \ q \ \text{erf}\left(\frac{q}{\sqrt{2}}\right) \exp\left(-\frac{\rho^2 (z-\beta)^2}{2 \delta_y^2 (1 + \rho^2 z^2)}\right) \right] }}
+with \eqn{ q = \frac{1 + \beta \rho^2 z}{\delta_y \sqrt{1 + \rho^2 z^2}} }
+and \eqn{ \text{erf}\left(\frac{q}{\sqrt{2}}\right) = \frac{2}{\sqrt{\pi}} \int_0^{\frac{q}{\sqrt{2}}}{\exp{(-t^2)}\ dt} }
+}
+\examples{
+# First example
+beta1 <- 0.15
+rho1 <- 5.75
+delta1 <- 0.22
+dz(0, bet = beta1, rho = rho1, delta = delta1)
+dz(0.5, bet = beta1, rho = rho1, delta = delta1)
+curve(dz(x, bet = beta1, rho = rho1, delta = delta1), from = -0.1, to = 0.7)
+
+# Second example
+beta2 <- 0.24
+rho2 <- 4.21
+delta2 <- 0.25
+dz(0, bet = beta2, rho = rho2, delta = delta2)
+dz(0.5, bet = beta2, rho = rho2, delta = delta2)
+curve(dz(x, bet = beta2, rho = rho2, delta = delta2), from = -0.1, to = 0.7)
+
+}
+\references{
+El Ghaziri, A., Bouhlel, N., Sapoukhina, N., Rousseau, D.,
+On the importance of non-Gaussianity in chlorophyll fluorescence imaging.
+Remote Sensing 15(2), 528 (2023).
+\doi{10.3390/rs15020528}
+
+Díaz-Francés, E., Rubio, F.J.,
+On the existence of a normal approximation to the distribution of the ratio of two independent normal random variables.
+Stat Papers 54, 309–323 (2013).
+\doi{10.1080/03610929808832115}
+}
+\author{
+Pierre Santagostini, Nizar Bouhlel
+}

man/kummerM.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/kummerM.R
+\name{kummerM}
+\alias{kummerM}
+\alias{lauricella}
+\title{Lauricella \eqn{D}-Hypergeometric Function}
+\usage{
+kummerM(a, b, z, eps = 1e-06)
+}
+\arguments{
+\item{a}{numeric.}
+
+\item{b}{numeric}
+
+\item{z}{numeric.}
+
+\item{eps}{numeric. Precision for the sum (default 1e-06).}
+}
+\value{
+A numeric value: the value of the Kummer's function,
+with two attributes \code{attr(, "epsilon")} (precision of the result) and \code{attr(, "k")} (number of iterations).
+}
+\description{
+Computes the Kummer's function, or confluent hypergeometric function.
+}
+\details{
+The Kummer's confluent hypergeometric function is given by:
+\deqn{\displaystyle{_1 F_1\left(a, b; z\right) = \sum_{n = 0}^{+\infty}{ \frac{ (a)_n }{ (b)_n } \frac{z^n}{n!} }}}
+
+where \eqn{(z)_p} is the Pochhammer symbol (see \code{\link{pochhammer}}).
+
+The \code{eps} argument gives the required precision for its computation.
+It is the \code{attr(, "epsilon")} attribute of the returned value.
+}
+\author{
+Pierre Santagostini, Nizar Bouhlel
+
+N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions.
+IEEE Signal Processing Letters, vol. 30, pp. 1672-1676, October 2023.
+\doi{10.1109/LSP.2023.3324594}
+}

man/lnpochhammer.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/lnpochhammer.R
+\name{lnpochhammer}
+\alias{lnpochhammer}
+\title{Logarithm of the Pochhammer Symbol}
+\usage{
+lnpochhammer(x, n)
+}
+\arguments{
+\item{x}{numeric.}
+
+\item{n}{positive integer.}
+}
+\value{
+Numeric value. The logarithm of the Pochhammer symbol.
+}
+\description{
+Computes the logarithm of the Pochhammer symbol.
+}
+\details{
+The Pochhammer symbol is given by:
+\deqn{ \displaystyle{ (x)_n = \frac{\Gamma(x+n)}{\Gamma(x)} = x (x+1) ... (x+n-1) } }
+So, if \eqn{n > 0}:
+\deqn{ \displaystyle{ log\left((x)_n\right) = log(x) + log(x+1) + ... + log(x+n-1) } }
+
+If \eqn{n = 0}, \eqn{\displaystyle{ log\left((x)_n\right) = log(1) = 0}}
+}
+\examples{
+lnpochhammer(2, 0)
+lnpochhammer(2, 1)
+lnpochhammer(2, 3)
+
+}
+\seealso{
+\code{\link{pochhammer}},
+\code{\link{lauricella}}
+}
+\author{
+Pierre Santagostini, Nizar Bouhlel
+}

man/pochhammer.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pochhammer.R
+\name{pochhammer}
+\alias{pochhammer}
+\title{Pochhammer Symbol}
+\usage{
+pochhammer(x, n)
+}
+\arguments{
+\item{x}{numeric.}
+
+\item{n}{positive integer.}
+}
+\value{
+Numeric value. The value of the Pochhammer symbol.
+}
+\description{
+Computes the Pochhammer symbol.
+}
+\details{
+The Pochhammer symbol is given by:
+\deqn{ \displaystyle{ (x)_n = \frac{\Gamma(x+n)}{\Gamma(x)} = x (x+1) ... (x+n-1) } }
+}
+\examples{
+pochhammer(2, 0)
+pochhammer(2, 1)
+pochhammer(2, 3)
+
+}
+\seealso{
+\code{\link{lauricella}}
+}
+\author{
+Pierre Santagostini, Nizar Bouhlel
+}