| Title: | Density, Cumulative and Quantile Functions of Quadratic Forms |
|---|---|
| Description: | The computation of quadratic form (QF) distributions is often not trivial and it requires numerical routines. The package contains functions aimed at evaluating the exact distribution of quadratic forms (QFs) and ratios of QFs. In particular, we propose to evaluate density, quantile and distribution functions of positive definite QFs and ratio of independent positive QFs by means of an algorithm based on the numerical inversion of Mellin transforms. |
| Authors: | Aldo Gardini [aut, cre], Fedele Greco [aut], Carlo Trivisano [aut] |
| Maintainer: | Aldo Gardini <[email protected]> |
| License: | GPL-3 |
| Version: | 0.0.6 |
| Built: | 2024-10-02 06:17:32 UTC |
| Source: | CRAN |
The function computes the Mellin transform of a positive definite quadratic form producing a MellinQF object.
The output can be used to evaluate the density, cumulative and quantile functions of
the target quadratic form.
compute_MellinQF( lambdas, etas = rep(0, length(lambdas)), eps = 1e-06, rho = 1 - 1e-04, maxit_comp = 1e+05, eps_quant = 1e-06, maxit_quant = 10000, lambdas_tol = NULL )compute_MellinQF( lambdas, etas = rep(0, length(lambdas)), eps = 1e-06, rho = 1 - 1e-04, maxit_comp = 1e+05, eps_quant = 1e-06, maxit_quant = 10000, lambdas_tol = NULL )
lambdas |
vector of positive weights. |
etas |
vector of non-centrality parameters. Default all zeros (central chi square). |
eps |
required absolute error for density and cumulative functions. |
rho |
distribution total probability mass for which it is desired to keep the error |
maxit_comp |
maximum number of iterations. |
eps_quant |
required numerical error for quantile computation. |
maxit_quant |
maximum number of iterations before stopping the quantile computation. |
lambdas_tol |
maximum value admitted for the weight skewness. When it is not NULL (default), elements of lambdas such that the ratio max(lambdas)/lambdas is greater than the specified value are removed. |
The quadratic form having positive weights lambdas and non-centrality parameters etas is considered:
Its Mellin transform is computed by exploiting the density formulation by Ruben (1962).
The numerical error is controlled in order to provide the requested precision (eps) for the
interval of quantiles that contains the specified total probability rho.
The argument eps_quant controls the relative precision requested for the
computation of quantiles that determine the range in which the error eps is
guaranteed, whereas maxit_quant sets the maximum number of Newton-Raphson iterations of the algorithm.
The function returns an object of the class MellinQF that contains information on the Mellin transform
of a linear combination of positively weighted chi-square random variables. This information can be used in order to
evaluate the density, cumulative distribution and quantile functions.
An object of the class MellinQF has the following components:
range_q: the range of quantiles that contains the specified mass of probability rho in which it
is possible to compute density and CDF preserving the error level eps.
Mellin: a list containing the values of the Mellin transform (Mellin),
the corresponding evaluation points (z), the integration step delta and the lowest weight (lambda_min).
the inputs rho, lambdas, etas, eps needed for CDF, PDF and quantile function computation.
Ruben, Harold. "Probability content of regions under spherical normal distributions, IV: The distribution of homogeneous and non-homogeneous quadratic functions of normal variables." The Annals of Mathematical Statistics 33.2 (1962): 542-570.
The function print.MellinQF can be used to summarize the basic information on the Mellin transform.
The object can be used in the function dQF to compute the density function of the QF,
pQF for the CDF and qQF for the quantile function.
library(QF) # Definition of the QF lambdas_QF <- c(rep(7, 6),rep(3, 2)) etas_QF <- c(rep(6, 6), rep(2, 2)) # Computation Mellin transform eps <- 1e-7 rho <- 0.999 Mellin <- compute_MellinQF(lambdas_QF, etas_QF, eps = eps, rho = rho) print(Mellin)library(QF) # Definition of the QF lambdas_QF <- c(rep(7, 6),rep(3, 2)) etas_QF <- c(rep(6, 6), rep(2, 2)) # Computation Mellin transform eps <- 1e-7 rho <- 0.999 Mellin <- compute_MellinQF(lambdas_QF, etas_QF, eps = eps, rho = rho) print(Mellin)
The function computes the Mellin transform of the ratio of independent and positive definite quadratic forms producing a MellinQF_ratio object.
The output can be used to evaluate the density, cumulative and quantile functions of the target quadratic form.
compute_MellinQF_ratio( lambdas_num, lambdas_den, etas_num = rep(0, length(lambdas_num)), etas_den = rep(0, length(lambdas_den)), eps = 1e-06, rho = 1 - 1e-04, maxit_comp = 1e+05, eps_quant = 1e-06, maxit_quant = 10000, lambdas_tol = NULL )compute_MellinQF_ratio( lambdas_num, lambdas_den, etas_num = rep(0, length(lambdas_num)), etas_den = rep(0, length(lambdas_den)), eps = 1e-06, rho = 1 - 1e-04, maxit_comp = 1e+05, eps_quant = 1e-06, maxit_quant = 10000, lambdas_tol = NULL )
lambdas_num |
vector of positive weights for the numerator. |
lambdas_den |
vector of positive weights for the denominator. |
etas_num |
vector of non-centrality parameters for the numerator. Default all zeros (central chi square). |
etas_den |
vector of non-centrality parameters for the denominator. Default all zeros (central chi square). |
eps |
required absolute error for density and cumulative functions. |
rho |
distribution total probability mass for which it is desired to keep the error |
maxit_comp |
maximum number of iterations. |
eps_quant |
required numerical error for quantile computation. |
maxit_quant |
maximum number of iterations before stopping the quantile computation. |
lambdas_tol |
maximum value admitted for the weight skewness for both the numerator and the denominator. When it is not NULL (default), elements of lambdas such that the ratio max(lambdas)/lambdas is greater than the specified value are removed. |
The Mellin transform of the ratio of two independent quadratic forms having positive weights lambdas_num and lambdas_den
and non-centrality parameters etas_num and etas_den is computed
exploiting the density formulation by Ruben (1962). The numerical error is controlled in order to provide the requested precision (eps) for the
interval of quantiles that contains the specified total probability rho.
The argument eps_quant controls the relative precision requested for the computation of quantiles that determine the range in which the error eps is
guaranteed, whereas maxit_quant sets the maximum number of Newton-Raphson iterations of the algorithm.
The function returns an object of the class MellinQF_ratio that contains information on the Mellin transform
of the ratio of two linear combinations of positively weighted chi-square random variables. This information can be used in order to
evaluate the density, cumulative distribution and quantile functions of the desired quadratic form.
An object of the class MellinQF_ratio has the following components:
range_q: the range of quantiles that contains the specified mass of probability rho in which it
is possible to compute density and CDF preserving the error level eps.
Mellin: a list containing the values of the Mellin transform (Mellin),
the corresponding evaluation points (z), the integration step delta and the lowest numerator weight (lambda_min).
the inputs rho, lambdas_num, lambdas_den, etas_num, etas_den, eps needed for CDF, PDF and quantile function computation.
Ruben, Harold. "Probability content of regions under spherical normal distributions, IV: The distribution of homogeneous and non-homogeneous quadratic functions of normal variables." The Annals of Mathematical Statistics 33.2 (1962): 542-570.
The function print.MellinQF_ratio can be used to summarize the basic information on the Mellin transform.
The object can be used in the function dQF_ratio to compute the density function of the QFs ratio,
pQF_ratio for the CDF and qQF_ratio for the quantile function.
library(QF) # Definition of the QFs lambdas_QF_num <- c(rep(7, 6),rep(3, 2)) etas_QF_num <- c(rep(6, 6), rep(2, 2)) lambdas_QF_den <- c(0.6, 0.3, 0.1) # Computation Mellin transform eps <- 1e-7 rho <- 0.999 Mellin_ratio <- compute_MellinQF_ratio(lambdas_QF_num, lambdas_QF_den, etas_QF_num, eps = eps, rho = rho) print(Mellin_ratio)library(QF) # Definition of the QFs lambdas_QF_num <- c(rep(7, 6),rep(3, 2)) etas_QF_num <- c(rep(6, 6), rep(2, 2)) lambdas_QF_den <- c(0.6, 0.3, 0.1) # Computation Mellin transform eps <- 1e-7 rho <- 0.999 Mellin_ratio <- compute_MellinQF_ratio(lambdas_QF_num, lambdas_QF_den, etas_QF_num, eps = eps, rho = rho) print(Mellin_ratio)
The function computes the CDF of the ratio of two dependent and possibly indefinite quadratic forms.
pQF_depratio( q = NULL, lambdas = NULL, A = NULL, B = NULL, eps = 1e-06, maxit_comp = 1e+05, lambdas_tol = NULL )pQF_depratio( q = NULL, lambdas = NULL, A = NULL, B = NULL, eps = 1e-06, maxit_comp = 1e+05, lambdas_tol = NULL )
q |
quantile. |
lambdas |
vector of eigenvalues of the matrix (A-qB). |
A |
matrix of the numerator QF. If not specified but |
B |
matrix of the numerator QF. If not specified but |
eps |
requested absolute error. |
maxit_comp |
maximum number of iterations. |
lambdas_tol |
maximum value admitted for the weight skewness for both the numerator and the denominator. When it is not NULL (default), elements of lambdas such that the ratio max(lambdas)/lambdas is greater than the specified value are removed. |
The distribution function of the following ratio of dependent quadratic forms is computed:
where .
The transformation to the following indefinite quadratic form is exploited:
The following inputs can be provided:
vector lambdas that contains the eigenvalues of the matrix . Input q is ignored.
matrix A and/or matrix B: in these cases q is required to be not null and an eventual missing specification of one matrix make it equal to the identity.
The values of the CDF at quantiles q.
MellinQF
It allows to visualize useful information on the MellinQF object.
## S3 method for class 'MellinQF' print(x, digits = 3L, ...)## S3 method for class 'MellinQF' print(x, digits = 3L, ...)
x |
|
digits |
number of digits to display in the output. |
... |
potentially more arguments passed to methods. |
MellinQF_ratio
It allows to visualize useful information on the MellinQF_ratio object.
## S3 method for class 'MellinQF_ratio' print(x, digits = 3L, ...)## S3 method for class 'MellinQF_ratio' print(x, digits = 3L, ...)
x |
|
digits |
number of digits to display in the output. |
... |
potentially more arguments passed to methods. |
Density function, distribution function, quantile function and random generator for positive definite QFs.
dQF(x, obj) pQF(q, obj) qQF(p, obj, eps_quant = 1e-06, maxit_quant = 10000) rQF(n, lambdas, etas = rep(0, length(lambdas)))dQF(x, obj) pQF(q, obj) qQF(p, obj, eps_quant = 1e-06, maxit_quant = 10000) rQF(n, lambdas, etas = rep(0, length(lambdas)))
x, q
|
vector of quantiles. |
obj |
|
p |
vector of probabilities. |
eps_quant |
relative error for quantiles. |
maxit_quant |
maximum number of Newton-Raphson iterations allowed to compute quantiles. |
n |
number of observations. |
lambdas |
vector of positive weights. |
etas |
vector of non-centrality parameters. Default all zeros. |
The quadratic form CDF and PDF are evaluated by numerical inversion of the Mellin transform.
The absolute error specified in compute_MellinQF is guaranteed for values of q and x inside the range_q.
If the quantile is outside range_q, computations are carried out, but a warning is sent.
The function uses the Newton-Raphson algorithm to compute the QF quantiles related to probabilities p.
dQF provides the values of the density function at a quantile x.
pQF provides the cumulative distribution function at a quantile q.
qQF provides the quantile corresponding to a probability level p.
rQF provides a sample of n independent realizations from the QF.
See compute_MellinQF for details on the Mellin computation.
library(QF) # Definition of the QF lambdas_QF <- c(rep(7, 6),rep(3, 2)) etas_QF <- c(rep(6, 6), rep(2, 2)) # Computation Mellin transform eps <- 1e-7 rho <- 0.999 Mellin <- compute_MellinQF(lambdas_QF, etas_QF, eps = eps, rho = rho) xs <- seq(Mellin$range_q[1], Mellin$range_q[2], l = 100) # PDF ds <- dQF(xs, Mellin) plot(xs, ds, type="l") # CDF ps <- pQF(xs, Mellin) plot(xs, ps, type="l") # Quantile qs <- qQF(ps, Mellin) plot(ps, qs, type="l") #Comparison computed quantiles vs real quantiles plot((qs - xs) / xs, type = "l")library(QF) # Definition of the QF lambdas_QF <- c(rep(7, 6),rep(3, 2)) etas_QF <- c(rep(6, 6), rep(2, 2)) # Computation Mellin transform eps <- 1e-7 rho <- 0.999 Mellin <- compute_MellinQF(lambdas_QF, etas_QF, eps = eps, rho = rho) xs <- seq(Mellin$range_q[1], Mellin$range_q[2], l = 100) # PDF ds <- dQF(xs, Mellin) plot(xs, ds, type="l") # CDF ps <- pQF(xs, Mellin) plot(xs, ps, type="l") # Quantile qs <- qQF(ps, Mellin) plot(ps, qs, type="l") #Comparison computed quantiles vs real quantiles plot((qs - xs) / xs, type = "l")
Density function, distribution function, quantile function and random generator for the ratio of positive definite QFs.
dQF_ratio(x, obj) pQF_ratio(q, obj) qQF_ratio(p, obj, eps_quant = 1e-06, maxit_quant = 10000) rQF_ratio( n, lambdas_num, lambdas_den, etas_num = rep(0, length(lambdas_num)), etas_den = rep(0, length(lambdas_den)) )dQF_ratio(x, obj) pQF_ratio(q, obj) qQF_ratio(p, obj, eps_quant = 1e-06, maxit_quant = 10000) rQF_ratio( n, lambdas_num, lambdas_den, etas_num = rep(0, length(lambdas_num)), etas_den = rep(0, length(lambdas_den)) )
x, q
|
vector of quantiles. |
obj |
|
p |
vector of probabilities. |
eps_quant |
relative error for quantiles. |
maxit_quant |
maximum number of Newton-Raphson iterations allowed to compute quantiles. |
n |
number of observations. |
lambdas_num |
vector of positive weights of the QF at the numerator. |
lambdas_den |
vector of positive weights of the QF at the denominator |
etas_num |
vector of non-centrality parameters of the QF at the numerator. Default all zeros. |
etas_den |
vector of non-centrality parameters of the QF at the denominator Default all zeros. |
The CDF and PDF of the ratio of positive QFs are evaluated by numerical inversion of the Mellin transform.
The absolute error specified in compute_MellinQF_ratio is guaranteed for values of q and x inside range_q.
If the quantile is outside range_q, computations are carried out, but a warning is sent.
The function uses the Newton-Raphson algorithm to compute the ratio of QFs quantiles related to probabilities p.
dQF_ratio provides the values of the density function at a quantile x.
pQF_ratio provides the cumulative distribution function at a quantile q.
qQF_ratio provides the quantile corresponding to a probability level p.
rQF_ratio provides a sample of n independent realizations the QFs ratio.
See compute_MellinQF_ratio for details on the Mellin computation.
lambdas_QF_num <- c(rep(7, 6),rep(3, 2)) etas_QF_num <- c(rep(6, 6), rep(2, 2)) lambdas_QF_den <- c(0.6, 0.3, 0.1) # Computation Mellin transform eps <- 1e-7 rho <- 0.999 Mellin_ratio <- compute_MellinQF_ratio(lambdas_QF_num, lambdas_QF_den, etas_QF_num, eps = eps, rho = rho) xs <- seq(Mellin_ratio$range_q[1], Mellin_ratio$range_q[2], l = 100) # PDF ds <- dQF_ratio(xs, Mellin_ratio) plot(xs, ds, type="l") # CDF ps <- pQF_ratio(xs, Mellin_ratio) plot(xs, ps, type="l") # Quantile qs <- qQF_ratio(ps, Mellin_ratio) plot(ps, qs, type="l") #Comparison computed quantiles vs real quantiles plot((qs - xs) / xs, type = "l")lambdas_QF_num <- c(rep(7, 6),rep(3, 2)) etas_QF_num <- c(rep(6, 6), rep(2, 2)) lambdas_QF_den <- c(0.6, 0.3, 0.1) # Computation Mellin transform eps <- 1e-7 rho <- 0.999 Mellin_ratio <- compute_MellinQF_ratio(lambdas_QF_num, lambdas_QF_den, etas_QF_num, eps = eps, rho = rho) xs <- seq(Mellin_ratio$range_q[1], Mellin_ratio$range_q[2], l = 100) # PDF ds <- dQF_ratio(xs, Mellin_ratio) plot(xs, ds, type="l") # CDF ps <- pQF_ratio(xs, Mellin_ratio) plot(xs, ps, type="l") # Quantile qs <- qQF_ratio(ps, Mellin_ratio) plot(ps, qs, type="l") #Comparison computed quantiles vs real quantiles plot((qs - xs) / xs, type = "l")