Estimation of the parameters of a ratio \(\displaystyle{Z = \frac{X}{Y}}\),
\(X\) and \(Y\) being two random variables distributed
according to Gaussian distributions,
using the EM (estimation-maximization) algorithm or variational inference.
Depending on the estimation method, the estparnormatio function calls
estparEM (EM algorithm) or estparVB (variational Bayes).
estparnormratio(z, method = c("EM", "VB"), indep = TRUE, eps = 1e-06,
na.rm = FALSE, display = FALSE, graph = FALSE,
xlim = NULL, ylim = NULL, mux0 = 1, sigmax0 = 1,
alphax0 = NULL, betax0 = NULL, muy0 = 1, sigmay0 = 1,
alphay0 = NULL, betay0 = NULL,
cov0 = 0)estparEM(z, indep = TRUE, eps = 1e-06, na.rm = FALSE,
display = FALSE, graph = FALSE, xlim = NULL, ylim = NULL,
mux0 = 1, sigmax0 = 1, muy0 = 1, sigmay0 = 1, cov0 = 0)
estparVB(z, eps = 1e-06, na.rm = FALSE, display = FALSE,
graph = FALSE, xlim = NULL, ylim = NULL,
mux0 = 1, sigmax0 = 1, alphax0 = 1, betax0 = 1,
muy0 = 1, sigmay0 = 1, alphay0 = 1, betay0 = 1)
A list of 4 elements beta, rho, delta, r:
the estimated parameters of the \(Z\) distribution
\(\hat{\beta}\), \(\hat{\rho}\), \(\hat{\delta}_y\) and \(\hat{r}\),
with three attributes attr(, "epsilon") (precision of the result),
attr(, "k") (number of iterations) and attr(, "method")
(estimation method).
If indep=FALSE, \(r\) is not estimated, it is set to 0.
numeric.
the method used to estimate the parameters of the distribution.
It can be "EM" (expectation-maximization) or "VB" (Variational Bayes).
logical. If indep=TRUE (default) \(X\) and \(Y\)
are two independent Gaussian variables and the parameters \(\beta\),
\(\rho\) and \(delta_y\) parameters of \(Z=\frac{X}{Y}\) are estimated.
If indep=FALSE, \(X\) and \(Y\) can be correlated, and
the correlation coefficient \(r\) is also estimated.
numeric. Precision for the estimation of the parameters (see Details).
a logical evaluating to TRUE or FALSE indicating
whether NA values should be stripped before the computation proceeds.
logical. When TRUE the successive values of the parameters
are printed.
logical. When TRUE the successive values of the parameters
are plotted.
if graph is TRUE, the x and y limits
of the plot. Default: xlim = c(0, 1000) and ylim depend on the initial
values of the parameters: ylim = 10*c(0, max(theta)).
initial values of the means and
standard deviations of the \(X\) and \(Y\) variables. Default:
mux0 = 1, sigmax0 = 1, muy0 = 1, sigmay0 = 1.
initial values for the variational
Bayes method. Omitted if method="EM".
If method="VB", if omitted, they are set to 1.
initial value of the covariance of \(X\) and \(Y\).
If indep is FALSE, cov0 must be different from 0.
Pierre Santagostini, Angélina El Ghaziri, Nizar Bouhlel
Let a random variable: \(\displaystyle{Z = \frac{X}{Y}}\),
\(X\) and \(Y\) being normally distributed: \(X \sim N(\mu_x, \sigma_x)\) and \(Y \sim N(\mu_y, \sigma_y)\).
The probability density of \(Z\) is: $$\displaystyle{ f_Z(z; \beta, \rho, \delta_y, r) = \frac{\rho \sqrt{1 - r^2}}{\pi (\rho^2 z^2 - 2r \rho z +1)} \ \exp{\left(-\frac{\rho^2 \beta^2 - 2r \beta \rho + 1}{2(1 - r^2) \delta_y^2}\right)} \ {}_1 F_1\left( 1, \frac{1}{2}; \frac{1}{2(1 - r^2) \delta_y^2} \frac{(\beta \rho^2 z - r\rho(z + \beta) + 1)^2}{\rho^2 z^2 - 2r \rho z + 1} \right) }$$
with: \(\displaystyle{\beta = \frac{\mu_x}{\mu_y}}\), \(\displaystyle{\rho = \frac{\sigma_y}{\sigma_x}}\), \(\displaystyle{\delta_y = \frac{\sigma_y}{\mu_y}}\), \(r = Cor(X, Y) = \frac{Cov(X, Y)}{\sigma_X \sigma_Y}\).
and \(_1 F_1\left(a, b; x\right)\) is the confluent hypergeometric function (Kummer's function): $$\displaystyle{ _1 F_1\left(a, b; x\right) = \sum_{n = 0}^{+\infty}{ \frac{ (a)_n }{ (b)_n } \frac{x^n}{n!} } }$$
If \(X\) and \(Y\) are independent (\(r = 0\)), the probability density is: $$\displaystyle{ f_Z(z; \beta, \rho, \delta_y, r = 0) = f_Z(z; \beta, \rho, \delta_y) = \frac{\rho}{\pi (1 + \rho^2 z^2)} \ \exp{\left(-\frac{\rho^2 \beta^2 + 1}{2\delta_y^2}\right)} \ {}_1 F_1\left( 1, \frac{1}{2}; \frac{1}{2 \delta_y^2} \frac{(1 + \beta \rho^2 z)^2}{1 + \rho^2 z^2} \right) }$$
If method = "EM", the means and standard deviations
\(\mu_x\), \(\sigma_x\), \(\mu_y\) and \(\sigma_y\)
and the correlation coefficient \(r\)
are estimated with the EM algorithm.
When \(r = 0\), \(\mu_x\), \(\sigma_x\), \(\mu_y\) and \(\sigma_y\) are estimated using the algorithm presented in El Ghaziri et al.
If method = "VB", they are estimated with the variational Bayes method
as presented in Bouhlel et al. For now, this method is available only
for the case when \(X\) and \(Y\) are independent, i.e. \(r = 0\).
Then the parameters \(\beta\), \(\rho\), \(\delta_y\) of the \(Z\) distribution are computed from these means and standard deviations.
The estimation of \(\mu_x\), \(\sigma_x\), \(\mu_y\) and \(\sigma_y\)
uses an iterative algorithm.
The precision for their estimation is given by the eps parameter.
The computation uses the kummer function.
If there are ties in the z vector, it generates a warning,
as there should be no ties in data distributed among a continuous distribution.
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). tools:::Rd_expr_doi("10.3390/rs15020528")
Bouhlel, N., Mercier, F., El Ghaziri, A., Rousseau, D., Parameter Estimation of the Normal Ratio Distribution with Variational Inference. 2023 31st European Signal Processing Conference (EUSIPCO), Helsinki, Finland, 2023, pp. 1823-1827. tools:::Rd_expr_doi("10.23919/EUSIPCO58844.2023.10290111")
dnormratio(): probability density of a normal ratio.
rnormratio(): sample simulation.
# \donttest{
# First example: ratio of independent variables
beta1 <- 0.15
rho1 <- 5.75
delta1 <- 0.22
set.seed(1234)
z1 <- rnormratio(800, bet = beta1, rho = rho1, delta = delta1)
# With the EM method:
estparnormratio(z1, method = "EM", indep = TRUE)
# With the variational method:
estparnormratio(z1, method = "VB")
# Second example: ratio of correlated variables
beta2 <- 0.24
rho2 <- 4.21
delta2 <- 0.25
r2 <- 0.8
set.seed(1234)
z2 <- rnormratio(800, bet = beta2, rho = rho2, r = r2, delta = delta2)
# With the EM method:
estparnormratio(z2, method = "EM", indep = FALSE)
# }
Run the code above in your browser using DataLab