README.Rmd

---
output: github_document
bibliography: vignettes/bibliography.bib
link-citations: TRUE
csl: vignettes/cran_style.csl
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```
```{r setup, include = FALSE}
require(qfratio)
set.seed(64501)
```

# qfratio: R Package for Moments and Distributions of Ratios of Quadratic Forms

<!-- badges: start -->
[![R-CMD-check](https://github.com/watanabe-j/qfratio/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/watanabe-j/qfratio/actions/workflows/R-CMD-check.yaml)
<!-- badges: end -->


This package provides functions to evaluate moments of ratios
(and products) of quadratic forms in normal variables, specifically using
recursive algorithms developed by @BaoKan2013 and @HillierEtAl2014.
Generating functions for these moments are closely related to the
top-order zonal and invariant polynomials of matrix arguments.
It also provides some functions to evaluate distribution and density
functions of simple ratios of quadratic forms in normal variables
using several methods from @Imhof1961, @Hillier2001,
@Forchini2002[@Forchini2005], @ButlerPaolella2008, and @BrodaPaolella2009.

There exist a couple of `Matlab` programs developed by Raymond Kan
(available from <https://www-2.rotman.utoronto.ca/~kan/>) for evaluating the
moments, but this `R` package is an independent project (not a fork or
translation) and has different functionalities, including evaluation of
moments of multiple ratios of a particular form and scaling to avoid numerical
overflow.
This has originally been developed for a biological application,
specifically for evaluating average evolvability measures in
evolutionary quantitative genetics [@Watanabe2023cevo], but can be used
for a broader class of statistics.



## Installation

***WARNING*** Installation size of this package can be very large
(>100 MB on Linux and macOS; ~3 MB on Windows with a recent version (`>= 4.2`)
of `Rtools`), as it involves lots of `RcppEigen` functions.


### From CRAN (stable version)
```{r eval = FALSE}
install.packages("qfratio")
```


### From GitHub (development version)
```{r eval = FALSE}
## Install devtools first:
# install.packages("devtools")

## Recommended installation (pandoc required):
devtools::install_github("watanabe-j/qfratio", dependencies = TRUE, build_vignettes = TRUE)

## Minimal installation:
# devtools::install_github("watanabe-j/qfratio")
```

### Dependencies

    Imports: Rcpp, MASS, stats
    LinkingTo: Rcpp, RcppEigen
    Suggests: mvtnorm, CompQuadForm, graphics, testthat (>= 3.0.0),
              knitr, rmarkdown

If installing from source, you also need [`pandoc`](https://pandoc.org) for
correctly building the vignette.
For `pandoc < 2.11`, `pandoc-citeproc` is required as well.
(Never mind if you use `RStudio`, which appears to have them bundled.)


## Examples

This package has two major functionalities: evaluating moments and
distribution function of ratios of quadratic forms in normal variables.

### Moments

This functionality concerns evaluation of the following moments:
$\mathrm{E} \left( \left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p /
\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q \right)$ and
$\mathrm{E} \left( \left( \mathbf{x}^T \mathbf{A} \mathbf{x} \right)^p /
\left( \mathbf{x}^T \mathbf{B} \mathbf{x} \right)^q
\left( \mathbf{x}^T \mathbf{D} \mathbf{x} \right)^r \right)$,
where $\mathbf{x} \sim N_n \left(\boldsymbol{\mu}, \boldsymbol{\Sigma}\right)$.

These quantities are evaluated by `qfrm(A, B, p, q, ...)` and
`qfmrm(A, B, D, p, q, r, ...)`.
Because they are evaluated as partial sums of infinite series
[@Smith1989; @Smith1993; @BaoKan2013; @HillierEtAl2009; @HillierEtAl2014],
the evaluation results come with an error bound (where available), and
a `plot` method is defined for inspecting numerical convergence.
```{r, examples}
## Simple matrices
nv <- 4
A <- diag(1:nv)
B <- diag(sqrt(nv:1))

## Expectation of (x^T A x)^2 / (x^T x)^2 where x ~ N(0, I)
qfrm(A, p = 2)

## Compare with Monte Carlo mean
mean(rqfr(1000, A = A, p = 2))

## Expectation of (x^T A x)^1/2 / (x^T x)^1/2
(mom_A0.5 <- qfrm(A, p = 1/2))

## Monte Carlo mean
mean(rqfr(1000, A = A, p = 1/2))

plot(mom_A0.5)

## Expectation of (x^T x) / (x^T A^-1 x)
##   = "average conditional evolvability"
(avr_cevoA <- qfrm(diag(nv), solve(A)))

mean(rqfr(1000, A = diag(nv), B = solve(A), p = 1))
plot(avr_cevoA)

## Expectation of (x^T x)^2 / (x^T A x) (x^T A^-1 x)
##   = "average autonomy"
(avr_autoA <- qfmrm(diag(nv), A, solve(A), p = 2, q = 1, r = 1))

mean(rqfmr(1000, A = diag(nv), B = A, D = solve(A), p = 2, q = 1, r = 1))
plot(avr_autoA)

## Expectation of (x^T A B x) / ((x^T A^2 x) (x^T B^2 x))^1/2
##   = "average response correlation"
## whose Monte Carlo evaluation is called the "random skewers" analysis,
## while this is essentially an analytic solution (with slight truncation error)
(avr_rcorA <- qfmrm(crossprod(A, B), crossprod(A), crossprod(B),
                    p = 1, q = 1/2, r = 1/2))

mean(rqfmr(1000, A = crossprod(A, B), B = crossprod(A), D = crossprod(B),
           p = 1, q = 1/2, r = 1/2))
plot(avr_rcorA)


## More complex (but arbitrary) example
## Expectation of (x^T A x)^2 / (x^T B x)^3 where x ~ N(mu, Sigma)
mu <- 1:nv / nv
Sigma <- diag(runif(nv) * 3)
(mom_A2B3 <- qfrm(A, B, p = 2, q = 3, mu = mu, Sigma = Sigma,
                  m = 500, use_cpp = TRUE))
plot(mom_A2B3)
```

### Distributions

This functionality concerns evaluation of the (cumulative) distribution
function, probability density, and quantiles of
$\left( \mathbf{x}^T \mathbf{A} \mathbf{x} /
        \mathbf{x}^T \mathbf{B} \mathbf{x} \right) ^ p$,
where $\mathbf{x} \sim N_n \left(\boldsymbol{\mu}, \boldsymbol{\Sigma}\right)$.

These are implemented in `pqfr(quantile, A, B, p, ...)`,
`dqfr(quantile, A, B, p, ...)`, and `qqfr(probability, A, B, p, ...)`,
 whose usage mimics that of regular distribution-related functions.
```{r, examples_distr}
## Example parameters
nv <- 4
A <- diag(1:nv)
B <- diag(sqrt(nv:1))
mu <- 1:nv * 0.2
quantiles <- 0:nv + 0.5

## Distribution function and density of
## (x^T A x) / (x^T B x) where x ~ N(0, I)
pqfr(quantiles, A, B)
dqfr(quantiles, A, B)

## 95, 99, and 99.9 percentiles of the same
qqfr(c(0.05, 0.01, 0.001), A, B, lower.tail = FALSE)

## Comparing profiles
qseq <- seq.int(1 / sqrt(nv) - 0.2, nv + 0.2, length.out = 100)

## Generate p-value sequences for
## (x^T A x) / (x^T B x) where x ~ N(0, I) vs
## (x^T A x) / (x^T B x) where x ~ N(mu, I)
pseq_central <- pqfr(qseq, A, B)
pseq_noncent <- pqfr(qseq, A, B, mu = mu)

## Graphical comparison
plot(qseq, type = "n", xlim = c(1 / sqrt(nv), nv), ylim = c(0, 1),
     xlab = "q", ylab = "F(q)")
lines(qseq, pseq_central, col = "royalblue4", lty = 1)
lines(qseq, pseq_noncent, col = "tomato", lty = 2)
legend("topleft", legend = c("central", "noncentral"),
       col = c("royalblue4", "tomato"), lty = 1:2)

## Generate density sequences for
## (x^T A x) / (x^T B x) where x ~ N(0, I) vs
## (x^T A x) / (x^T B x) where x ~ N(mu, I)
dseq_central <- dqfr(qseq, A, B)
dseq_noncent <- dqfr(qseq, A, B, mu = mu)

## Graphical comparison
plot(qseq, type = "n", xlim = c(1 / sqrt(nv), nv), ylim = c(0, 0.7),
     xlab = "q", ylab = "f(q)")
lines(qseq, dseq_central, col = "royalblue4", lty = 1)
lines(qseq, dseq_noncent, col = "tomato", lty = 2)
legend("topright", legend = c("central", "noncentral"),
       col = c("royalblue4", "tomato"), lty = 1:2)
```


## Copyright notice

This package bundles selected `C` codes and part of `config.ac` from the
[GNU Scientific Library](https://www.gnu.org/software/gsl/),
whose copyright belongs to the original authors.
See `DESCRIPTION` and individual code files in `src/gsl` for details.
The redistribution complies with the GNU General Public License version 3.


## References