Package 'scorematchingad'

Description

Hyvärinen's score matching (Hyvärinen, 2005) <https://jmlr.org/papers/v6/hyvarinen05a.html> is a useful estimation technique when the normalising constant for a probability distribution is difficult to compute. This package implements score matching estimators using automatic differentiation in the 'CppAD' library <https://github.com/coin-or/CppAD> and is designed for quickly implementing score matching estimators for new models. Also available is general robustification (Windham, 1995) <https://www.jstor.org/stable/2346159>. Already in the package are estimators for directional distributions (Mardia, Kent and Laha, 2016) <doi:10.48550/arXiv.1604.08470> and the flexible Polynomially-Tilted Pairwise Interaction model for compositional data. The latter estimators perform well when there are zeros in the compositions (Scealy and Wood, 2023) <doi:10.1080/01621459.2021.2016422>, even many zeros (Scealy, Hingee, Kent, and Wood, 2024) <doi:10.1007/s11222-024-10412-w>.

Details

This package's main features are

• A general capacity to implement score matching estimators that use algorithmic differentiation to avoid tedious manual algebra. The package uses CppAD and Eigen to differentiate model densities and compute the score matching discrepancy function (see scorematchingtheory). The score matching discrepancy is usually minimised by solving a quadratic equation, but a method for solving numerically (through optimx::Rcgmin()) is also included. On Linux platforms using the gcc compiler new models can be fitted with the help of customll(), in a similar fashion to models in the TMB package. New manifolds or new transforms require small alterations to the source code of this package.

• Score matching estimators for the Polynomially-Tilted Pairwise Interaction (PPI) model (Scealy and Wood 2023; Scealy et al. 2024). See function ppi().

• Score matching and hybrid score matching estimators for von Mises Fisher, Bingham and Fisher-Bingham directional distributions (Mardia et al. 2016). See vMF(), Bingham() and FB().

• Implementation of a modification of Windham's robustifying method (Windham 1995) for many exponential family distributions. See Windham(). For some models the density approaches infinity at some locations, creating difficulties for the weights in Windham's original method (Scealy et al. 2024).

Acknowledgements

Colleagues Andrew T. A. Wood and John T. Kent played important roles in developing the statistical ideas and theory for score matching estimation for the PPI model (Scealy et al. 2024).

We developed this package on Ngunnawal and Ngambri Country. We thank the Country for its influence.

Author(s)

Maintainer: Kassel Liam Hingee [email protected] (ORCID)

Authors:

• Janice Scealy (ORCID)

Other contributors:

• Bradley M. Bell [copyright holder]

References

Amaral GJA, Dryden IL, Wood ATA (2007). “Pivotal Bootstrap Methods for k-Sample Problems in Directional Statistics and Shape Analysis.” Journal of the American Statistical Association, 102(478), 695–707. 27639898, http://www.jstor.org/stable/27639898.

Bell B (2023). “CppAD: A Package for Differentiation of C++ Algorithms.” https://github.com/coin-or/CppAD.

Hyvärinen A (2005). “Estimation of Non-Normalized Statistical Models by Score Matching.” Journal of Machine Learning Research, 6(24), 695–709. https://jmlr.org/papers/v6/hyvarinen05a.html.

Hyvärinen A (2007). “Some extensions of score matching.” Computational Statistics & Data Analysis, 51(5), 2499–2512. doi:10.1016/j.csda.2006.09.003.

Liu S, Kanamori T, Williams DJ (2019). “Estimating Density Models with Truncation Boundaries using Score Matching.” doi:10.48550/arXiv.1910.03834.

Mardia K (2018). “A New Estimation Methodology for Standard Directional Distributions.” In 2018 21st International Conference on Information Fusion (FUSION), 724–729. doi:10.23919/ICIF.2018.8455640.

Mardia KV, Jupp PE (2000). Directional Statistics, Probability and Statistics. Wiley, Great Britain. ISBN 0-471-95333-4.

Mardia KV, Kent JT, Laha AK (2016). “Score matching estimators for directional distributions.” doi:10.48550/arXiv.1604.08470.

Martin I, Uh H, Supali T, Mitreva M, Houwing-Duistermaat JJ (2019). “The mixed model for the analysis of a repeated-measurement multivariate count data.” Statistics in Medicine, 38(12), 2248–2268. doi:10.1002/sim.8101.

Scealy JL, Hingee KL, Kent JT, Wood ATA (2024). “Robust score matching for compositional data.” Statistics and Computing, 34, 93. doi:10.1007/s11222-024-10412-w.

Scealy JL, Wood ATA (2023). “Score matching for compositional distributions.” Journal of the American Statistical Association, 118(543), 1811–1823. doi:10.1080/01621459.2021.2016422.

Windham MP (1995). “Robustifying Model Fitting.” Journal of the Royal Statistical Society. Series B (Methodological), 57(3), 599–609. 2346159, http://www.jstor.org/stable/2346159.

Wiria AE, Prasetyani MA, Hamid F, Wammes LJ, Lell B, Ariawan I, Uh HW, Wibowo H, Djuardi Y, Wahyuni S, Sutanto I, May L, Luty AJ, Verweij JJ, Sartono E, Yazdanbakhsh M, Supali T (2010). “Does treatment of intestinal helminth infections influence malaria?” BMC Infectious Diseases, 10, 77. doi:10.1186/1471-2334-10-77.

Yu S, Drton M, Shojaie A (2019). “Generalized Score Matching for Non-Negative Data.” Journal of Machine Learning Research, 20(76), 1–70. https://jmlr.org/papers/v20/18-278.html.

Yu S, Drton M, Shojaie A (2020). “Generalized Score Matching for General Domains.” doi:10.48550/arXiv.2009.11428.

See Also

Useful links:

• https://github.com/kasselhingee/scorematchingad

• Report bugs at https://github.com/kasselhingee/scorematchingad/issues

Description

This is a low level object useful for implementing score matching estimators. An R6 class for storing a 'pointer' to a CppAD tape in ⁠C++⁠ (also called an ADFun) and associated information. Currently tools for modifying this information are not available in this package, however tools for creating new ADFun objects from an existing ADFun are available. Typically an ADFun object will be created by buildsmdtape().

Introduction to CppAD Tapes

This package uses version 2024000.5 of the algorithmic differentiation library CppAD (Bell 2023) to build score matching estimators. Full help for CppAD can be found at https://cppad.readthedocs.io/.

Differentiation proceeds by taping the basic (atomic) operations performed on the independent variables and dynamic parameters. The atomic operations include multiplication, division, addition, sine, cosine, exponential and many more. Example values for the variables and parameters are used to conduct this taping, so care must be taken with any conditional (e.g. if-then) operations, and CppAD has a special tool for this called CondExp (short for ⁠conditional expressions⁠). The result of taping is an object of class ADFun in CppAD and is often called a tape. This ADFun object can be evaluated, differentiated, used for further taping (via CppAD's base2ad()), solving differential equations and more. The differentiation is with respect to the independent variables, however the dynamic parameters can be altered which allows for creating a new ADFun object where the dynamic parameters become independent variables (see tapeSwap()). For the purposes of score matching, there are also fixed parameters, which are the elements of the model's parameter vector that are given and not estimated.

Warning: multiple CPU

Each time a tape is evaluated the corresponding ⁠C++⁠ object is altered. Parallel use of the same ADFun object thus requires care and is not tested. For now I recommend creating a new ADFun object for each CPU.

Active bindings

Methods

Public methods

• ADFun$new()

Method new()

Create a new ADFun object from an external pointer to a CppAD ADFun object.

Usage

ADFun$new(
  ptr,
  name = NULL,
  xtape = vector("numeric"),
  dyntape = vector("numeric"),
  usertheta = rep(NA_real_, length(dyntape))
)

Arguments

Examples

tapes <- buildsmdtape(
  "sim", "sqrt", "sph",
  ll = "ppi",
  ytape =  rep(1/3, 3),
  usertheta = ppi_paramvec(p=3),
  bdryw = "minsq",
  acut = 0.01,
  verbose = FALSE)
tapes$smdtape$xtape
tapes$smdtape$dyntape
tapes$smdtape$name
tapes$smdtape$ptr

Description

Score matching estimators for the Bingham distribution's parameter matrix. Two methods are available: a full score matching method that estimates the parameter matrix directly and a hybrid method by Mardia et al. (2016) that uses score matching to estimate just the eigenvalues of the parameter matrix.

Usage

Bingham(Y, A = NULL, w = rep(1, nrow(Y)), method = "Mardia")

Arguments

Details

The Bingham distribution has a density proportional to

$\exp(z^T A z),$

where $A$ is a symmetric matrix and the trace (sum of the diagonals) of $A$ is zero for identifiability (p181, Mardia and Jupp 2000).

The full score matching method estimates all elements of $A$ directly except the final element of the diagonal, which is calculated from the sum of the other diagonal elements to ensure that the trace of $A$ is zero.

The method by Mardia et al. (2016) first calculates the maximum-likelihood estimate of the eigenvectors $G$ of $A$. The observations Y are then standardised to Y$G$. This standardisation corresponds to diagonalising $A$ where the eigenvalues of $A$ become the diagonal elements of the new $A$. The diagonal elements of the new $A$ are then estimated using score matching, with the final diagonal element calculated from the sum of the other elements. See Mardia et al. (2016) for details.

Value

A list of est, SE and info.

• est contains the estimated matrix A and a vector form, paramvec, of A (ordered according to c(diag(A)[1:(p-1)], A[upper.tri(A)]) ). For the Mardia method, the estimated eigenvalues of A (named evals) and eigenvectors of A (named G) are also returned.

• SE contains estimates of the standard errors if computed. See cppad_closed().

• info contains a variety of information about the model fitting procedure and results.

References

Mardia KV, Jupp PE (2000). Directional Statistics, Probability and Statistics. Wiley, Great Britain. ISBN 0-471-95333-4.

Mardia KV, Kent JT, Laha AK (2016). “Score matching estimators for directional distributions.” doi:10.48550/arXiv.1604.08470.

See Also

Other directional model estimators: FB(), vMF_robust(), vMF()

Examples

p <- 4
A <- rsymmetricmatrix(p)
A[p,p] <- -sum(diag(A)[1:(p-1)]) #to satisfy the trace = 0 constraint
if (requireNamespace("simdd")){
  Y <- simdd::rBingham(100, A)
  Bingham(Y, method = "Mardia")
}

Description

For a parametric model family, the function buildsmdtape() generates CppAD tapes (called ADFuns) for the improper log-likelihood (without normalising constant) of the family and the score matching discrepancy function $A(z) + B(z) + C(z)$ (defined in scorematchingtheory). Three steps are performed by buildsmdtape(): first an object that specifies the manifold and any transformation to another manifold is created; then a tape of the log-likelihood (without normalising constant) is created; finally a tape of $A(z) + B(z) + C(z)$ is created.

Usage

buildsmdtape(
  start,
  tran = "identity",
  end = start,
  ll,
  ytape,
  usertheta,
  bdryw = "ones",
  acut = 1,
  thetatape_creator = function(n) {     seq(length.out = n) },
  verbose = FALSE
)

Arguments

Details

The improper log-likelihood (without normalising constant) must be implemented in ⁠C++⁠ and is selected by name. Similarly the transforms of the manifold must be implemented in ⁠C++⁠ and selected by name.

When using, CppAD one first creates tapes of functions. These tapes can then be used for evaluating the function and its derivatives, and generating further tapes through argument swapping, differentiation and composition. The taping relies on specifying typical argument values for the functions (see Introduction to CppAD Tapes below). Tapes can have both independent variables and dynamic parameters. The differentiation with CppAD occurs with respect to the independent variables. Tapes of tapes are possible, including tapes that swap the independent and dynamic variables - this is how this package differentiates with respect to a dynamic variables (see tapeSwap()).

To build a tape for the score matching discrepancy function, the package first tapes the map from a point $z$ on the end manifold to the value of the improper log-likelihood, where the independent variable is the $z$, the dynamic parameter is a vector of the parameters to estimate, and the remaining model parameters are fixed and not estimated. This tape is then used to generate a tape for the score matching discrepancy function where the parameters to estimate are the independent variable.

Only some combinations of start, tran and end are available because tran must map between start and end. These combinations of start-tran-end are currently available:

• sim-sqrt-sph

• sim-identity-sim

• sim-alr-Euc

• sim-clr-Hn111

• sph-identity-sph

• Euc-identity-Euc

Currently available improper log-likelihood functions are:

• dirichlet

• ppi

• vMF

• Bingham

• FB

Value

A list of:

• an ADFun object containing a tape of an improper likelihood with $z$ on the end manifold as the independent variable

• an ADFun object containing a tape of the score matching discrepancy function with the non-fixed parameters as the independent variable, and the measurements on the end manifold as the dynamic parameter.

• some information about the tapes

Introduction to CppAD Tapes

This package uses version 2024000.5 of the algorithmic differentiation library CppAD (Bell 2023) to build score matching estimators. Full help for CppAD can be found at https://cppad.readthedocs.io/.

Differentiation proceeds by taping the basic (atomic) operations performed on the independent variables and dynamic parameters. The atomic operations include multiplication, division, addition, sine, cosine, exponential and many more. Example values for the variables and parameters are used to conduct this taping, so care must be taken with any conditional (e.g. if-then) operations, and CppAD has a special tool for this called CondExp (short for ⁠conditional expressions⁠). The result of taping is an object of class ADFun in CppAD and is often called a tape. This ADFun object can be evaluated, differentiated, used for further taping (via CppAD's base2ad()), solving differential equations and more. The differentiation is with respect to the independent variables, however the dynamic parameters can be altered which allows for creating a new ADFun object where the dynamic parameters become independent variables (see tapeSwap()). For the purposes of score matching, there are also fixed parameters, which are the elements of the model's parameter vector that are given and not estimated.

Warning: multiple CPU

Each time a tape is evaluated the corresponding ⁠C++⁠ object is altered. Parallel use of the same ADFun object thus requires care and is not tested. For now I recommend creating a new ADFun object for each CPU.

Warning

There is no checking of the inputs ytape and usertheta.

References

Bell B (2023). “CppAD: A Package for Differentiation of C++ Algorithms.” https://github.com/coin-or/CppAD.

See Also

Other tape builders: moretapebuilders

Examples

p <- 3
u <- rep(1/sqrt(p), p)
ltheta <- p #length of vMF parameter vector
intheta <- rep(NA, length.out = ltheta)
tapes <- buildsmdtape("sph", "identity", "sph", "vMF",
              ytape = u,
              usertheta = intheta,
              "ones", verbose = FALSE
              )
evaltape(tapes$lltape, u, runif(n = ltheta))
evaltape(tapes$smdtape, runif(n = ltheta), u)

u <- rep(1/3, 3)
tapes <- buildsmdtape("sim", "sqrt", "sph", "ppi",
              ytape = u,
              usertheta = ppi_paramvec(p = 3),
              bdryw = "minsq", acut = 0.01,
              verbose = FALSE
              )
evaltape(tapes$lltape, u, rppi_egmodel(1)$theta)
evaltape(tapes$smdtape, rppi_egmodel(1)$theta, u)

Description

For a tape (i.e. an ADFun object) of a quadratic-form score matching discrepancy function, calculates the vector of parameters such that the gradient of the score matching discrepancy is zero. Also estimates standard errors and covariance. Many score matching estimators have an objective function that has a quadratic form (see scorematchingtheory).

Usage

cppad_closed(
  smdtape,
  Y,
  Yapproxcentres = NA * Y,
  w = rep(1, nrow(Y)),
  approxorder = 10
)

Arguments

Details

When the score matching discrepancy function is of quadratic form, then the gradient of the score matching discrepancy is zero at $H^{-1}b$, where $H$ is the average of the Hessian of the score matching discrepancy function evaluated at each measurement and $b$ is the average of the gradient offset (see quadratictape_parts()) evaluated at each measurement. Both the Hessian and the gradient offset are constant with respect to the model parameters for quadratic-form score matching discrepancy functions.

Standard errors use the Godambe information matrix (aka sandwich method) and are only computed when the weights are constant. The estimate of the negative of the sensitivity matrix $-G$ is the average of the Hessian of smdtape evaluated at each observation in Y. The estimate of the variability matrix $J$ is the sample covariance (denominator of $n-1$) of the gradient of smdtape evaluated at each of the observations in Y for the estimated $\theta$. The estimated variance of the estimator is then as $G^{-1}JG^{-1}/n,$ where n is the number of observations.

Taylor approximation is available because boundary weight functions and transformations of the measure in Hyvärinen divergence can remove singularities in the model log-likelihood, however evaluation at these singularities may still involve computing intermediate values that are unbounded. If the singularity is ultimately removed, then Taylor approximation from a nearby location will give a very accurate evaluation at the removed singularity.

See Also

Other generic score matching tools: Windham(), cppad_search()

Examples

smdtape <- buildsmdtape("sim", "sqrt", "sph", "ppi",
              ytape = rep(1/3, 3),
              usertheta = ppi_paramvec(p = 3),
              bdryw = "minsq", acut = 0.01,
              verbose = FALSE
              )$smdtape
Y <- rppi_egmodel(100)
cppad_closed(smdtape, Y$sample)

Description

Uses conjugate gradient descent to search for a vector of parameters such that gradient of the score matching discrepancy is within tolerance of zero. Also estimates standard errors and covariance.

Usage

cppad_search(
  smdtape,
  theta,
  Y,
  Yapproxcentres = NA * Y,
  w = rep(1, nrow(Y)),
  approxorder = 10,
  control = list(tol = 1e-15, checkgrad = TRUE)
)

Arguments

Details

The score matching discrepancy function and gradient of the score matching function are passed to optimx::Rcgmin(). The call to optimx::Rcgmin() uses the sum of observations (as opposed to the mean) to reduce floating point inaccuracies. This has implications for the meaning of the control parameters passed to Rcgmin() (e.g. tol). The results are converted into averages so the use of sums can be ignored when not setting control parameters, or studying the behaviour of Rcgmin.

Standard errors use the Godambe information matrix (aka sandwich method) and are only computed when the weights are constant. The estimate of the sensitivity matrix $G$ is the negative of the average over the Hessian of smdtape evaluated at each observation in Y. The estimate of the variability matrix $J$ is then the sample covariance (denominator of $n-1$) of the gradient of smdtape evaluated at each of the observations in Y for the estimated $\theta$. The variance of the estimator is then estimated as $G^{-1}JG^{-1}/n,$ where n is the number of observations.

Taylor approximation is available because boundary weight functions and transformations of the measure in Hyvärinen divergence can remove singularities in the model log-likelihood, however evaluation at these singularities may still involve computing intermediate values that are unbounded. If the singularity is ultimately removed, then Taylor approximation from a nearby location will give a very accurate evaluation at the removed singularity.

Warning

There appears to be floating point issues with evaluation of the gradient leading to slow or no convergence in many situations. Tweaking the convergence tolerance tol may be appropriate. If the closed form solution exists please use cppad_closed() instead.

See Also

Other generic score matching tools: Windham(), cppad_closed()

Examples

smdtape <- buildsmdtape("sim", "sqrt", "sph", "ppi",
              ytape = rep(1/3, 3),
              usertheta = ppi_paramvec(p = 3),
              bdryw = "minsq", acut = 0.01,
              verbose = FALSE
              )$smdtape
Y <- rppi_egmodel(100)
cppad_search(smdtape, 0.9 * Y$theta, Y$sample)
sum((smvalues_wsum(smdtape, Y$sample, Y$theta)$grad/nrow(Y$sample))^2)

Description

Supply ⁠C++⁠ code to specify a custom log-likelihood, much like TMB::compile() is passed ⁠C++⁠ code that formulate models. For score matching the normalising constant of the log-likelihood can be omitted. Taping of the function currently only works with the gcc compiler, typically on Linux; taping works if customll_test() returns TRUE. Packages RcppEigen and RcppXPtrUtils must be installed.

Usage

customll(
  code,
  rebuild = FALSE,
  includes = character(),
  cacheDir = getOption("rcpp.cache.dir", tempdir()),
  showOutput = verbose,
  verbose = getOption("verbose")
)

customll_test()

Arguments

Details

The function uses RcppXPtrUtils::cppXPtr() and Rcpp::cppFunction(). It is good practice to check the returned object using evalll(). The first compilation in a session can be very slow.

Value

customll() returns an adloglikelood object (which is just an externalptr with attributes) for the compiled log-likelihood function. The returned object has an attribute fname.

customll_test() returns TRUE if taping works in your R session. Otherwise it will return FALSE and generate warnings.

Code Argument

code must be ⁠C++⁠ that uses only CppAD and Eigen, which makes it very similar to the requirements of the input to TMB::compile() (which also uses CppAD and Eigen).

The start of code should always be "⁠a1type fname(const veca1 &x, const veca1 &theta){⁠" where fname is your chosen name of the log-likelihood function, x represents a point in the data space and theta is a vector of parameters for the log-likelihood. This specifies that the function will have two vector arguments (of type veca1) and will return a single numeric value (a1type).

The type a1type is a double with special ability for being taped by CppAD. The veca1 type is a vector of a1type elements, with the vector wrapping supplied by the Eigen C++ package (that is an Eigen matrix with 1 column and dynamic number of rows).

The body of the function must use operations from Eigen and/or CppAD, prefixed by ⁠Eigen::⁠ and ⁠CppAD::⁠ respectively. There are no easy instructions for writing these as it is genuine ⁠C++⁠ code, which can be very opaque to those unfamiliar with ⁠C++⁠. I've found the quick reference pages for for Eigen useful. Limited unary and binray operations are available directly from CppAD without Eigen. For the purposes of score matching the operations should all be smooth to create a smooth log-likelihood and the normalising constant may be omitted.

Examples

## Not run: customll_test()

myll <- customll("a1type dirichlet(const veca1 &u, const veca1 &beta) {
  size_t d  = u.size();
  a1type y(0.);  // initialize summation at 0
  for(size_t i = 0; i < d; i++)
  {   y   += beta[i] * log(u[i]);
  }
  return y;
}")
evalll(myll, rep(1/3, 3), rep(-0.5, 3))

tapes <- buildsmdtape("sim", "identity", "sim", 
 myll, rep(1/3, 3), rep(NA, 3), 
 bdryw="minsq", acut = 0.01)
evaltape(tapes$lltape, rep(1/3, 3), rep(-0.5, 3))
## End(Not run)

Description

Compute the natural logarithm of the improper density for the PPI model for the given matrix of measurements Y. Rows with negative values or with a sum that differs from 1 by more than 1E-15 are assigned a value of -Inf.

Usage

dppi(Y, ..., paramvec = NULL)

Arguments

Details

The value calculated by dppi is

$z_L^TA_Lz_L + b_L^Tz_L + \beta^T \log(z),$

where $z$ is the multivariate observation (i.e. a row of Y), and $z_L$ omits the final element of $z$.

See Also

Other PPI model tools: ppi_param_tools, ppi_robust(), ppi(), rppi()

Examples

m <- rppi_egmodel(10)
dppi(m$sample, paramvec = m$theta)

Description

Evaluates a custom log-likelihood function from customll() without taping. This is useful to check that the custom log-likelihood is behaving. To check a tape of the custom log-likelihood use buildsmdtape() then evaltape().

Usage

evalll(ll, u, theta)

Arguments

Value

The value of the log-likelihood at u with parameters theta.

Description

Evaluates a tape exactly or approximately for an array of provided variable values and dynamic parameter values. The function evaltape_wsum() computes the column-wise weighted sum of the result.

Usage

evaltape(tape, xmat, pmat, xcentres = NA * xmat, approxorder = 10)

evaltape_wsum(
  tape,
  xmat,
  pmat,
  w = NULL,
  xcentres = NA * xmat,
  approxorder = 10
)

Arguments

Details

Approximation is via Taylor approximation of the independent variable around the approximation centre provided in xcentres.

Value

A matrix, each row corresponding to the evaluation of the same row in xmat, pmat and xcentres.

See Also

Other tape evaluators: quadratictape_parts(), smvalues(), testquadratic()

Examples

u <- rep(1/3, 3)
tapes <- buildsmdtape("sim", "sqrt", "sph", "ppi",
              ytape = u,
              usertheta = ppi_paramvec(p = 3),
              bdryw = "minsq", acut = 0.01,
              verbose = FALSE
              )
evaltape(tapes$lltape, u, rppi_egmodel(1)$theta)
evaltape(tapes$smdtape, rppi_egmodel(1)$theta, u)
evaltape(tapes$lltape, rbind(c(0, 0, 1), c(0,0,1)), 
         rppi_egmodel(1)$theta, 
         xcentres = rbind(c(0.0005, 0.0005, 0.999), NA))

Description

Estimates parameters for the Fisher-Bingham distribution using score-matching.

Usage

FB(Y, km = NULL, A = NULL)

Arguments

Details

The density of the Fisher-Bingham distribution is proportional to

$\exp(z^TAz + \kappa\mu^Tz),$

where $A$ is a matrix as in the Bingham distribution, and $\kappa$ and $\mu$ are the concentration and mean direction, respectively, as in the von Mises-Fisher distribution.

Warning: Slow Convergence with Sample Size

Score matching estimates of all elements of $A$ and $\kappa\mu$ converge slowly with sample size. Even with a million simulated measurements, the gradient of the score matching discrepancy at the true parameters can have size (L2 Euclidean norm) more than 0.001, which is substantially non-zero.

See Also

Other directional model estimators: Bingham(), vMF_robust(), vMF()

Examples

p <- 3
A <- rsymmetricmatrix(p, -10, 10)
A[p,p] <- -sum(diag(A)[1:(p-1)]) #to satisfy the trace = 0 constraint
m <- runif(p, -10, 10)
m <- m / sqrt(sum(m^2))
if (requireNamespace("simdd")){
  Y <- simdd::rFisherBingham(1000, 2 * m, A)
  FB(Y)
}

Description

The microbiome data contains paired DNA samples from before treatment and 21 months after treatment for helminth infections (Martin et al. 2019). This data was analysed by Martin et al. (2019) and a further subset was studied by Scealy and Wood (2023). The data are from a study into the effect of helminth infections on the course of malaria infections (ImmunoSPIN-Malaria) in the Nangapanda subdistrict, Indonesia (Wiria et al. 2010). As part of the study, some participants were given 400mg of albendazole every three months for 1.5 years, remaining participants were given a placebo (Wiria et al. 2010).

Usage

microbiome

Format

A dataframe with 300 rows (two rows per individual) and 31 columns:

IndividualID

An integer uniquely specifying the individual.

Year

The collection year for the sample. 2008 for before treatment. 2010 for after treatment.

Sex

1 if female, 0 otherwise.

Treatment

TRUE if individual given 400mg of albendazole every three months for 1.5 years, FALSE otherwise.

Age

Age at first sample.

ct_Al

A Helminth measurement: The qPCR cycle threshold (CT) for Ascaris lumbricoides (large roundworm). Ascaris lumbricoides can be considered present if the value is 30 or less.

ct_Na

A Helminth measurement: The qPCR cycle threshold (CT) for Necator americanus (a hookworm). Necator americanus can be considered present if the value is 30 or less.

ct_Ad

A Helminth measurement: The qPCR cycle threshold (CT) for Ancylostoma duodenale (a hookworm). Ancylostoma duodenale can be considered present if the value is 30 or less.

micr_Tt

A Helminth measurement: The presence of Trichuris trichiura as determined by microscopy. A value of TRUE means Trichuris trichiura was detected.

Helminth

A Helminth measurement: If any of the above helminths were detected then TRUE, otherwise FALSE.

Remaining columns

Count prevalence of 18 bacterial phyla and 2 unclassified columns.

Details

The measurements in the data come from stool samples before and after treatment. Gut microbiome prevalence was measured using 16s rRNA 454 sequencing (Martin et al. 2019). Helminth infections were detected by PCR or microscopy (Martin et al. 2019).

The subset studied by Scealy and Wood (2023) contained only the measurements from before treatment, and only those individuals with a helminth infection. These measurements can be obtained by running

microbiome[(microbiome$Year == 2008) & microbiome$Helminth, ]

Two further individuals (IndividualID of 2079 and 2280) were deemed outliers by Scealy and Wood (2023).

Modifications from the Source

The microbiome data was created from the file S1_Table.xlsx hosted on Nematode.net at ⁠http://nematode.net/Data/environmental_interaction/S1_Table.xlsx⁠ using the below code.

microbiome <- readxl::read_excel("S1_Table.xlsx",
  range = "A3:AE303") #avoids the genus data, keeping - only phyla
metacolnames <- readxl::read_excel("S1_Table.xlsx",
  range = "A2:J2", 
  col_names = FALSE)
colnames(microbiome)[1:ncol(metacolnames)] <- metacolnames[1, ]
colnames(microbiome)[2] <- "Year"
microbiome[, 11] <- (microbiome$ct_Al <= 30) | (microbiome$ct_Na <= 30) |
  (microbiome$ct_Ad <= 30) | (microbiome$ct_St <= 30) |
  (microbiome$micr_Tt == 1)
colnames(microbiome)[11] <- "Helminth"
microbiome <- microbiome |>
  dplyr::mutate(across(c(1,2,3,12:31), as.integer)) |>
  dplyr::mutate(micr_Tt = as.logical(micr_Tt),
                Treatment = as.logical(Treatment)) |>
  dplyr::rename(IndividualID = `Individual ID`)
microbiome <- as.data.frame(microbiome)

Source

⁠http://nematode.net/Data/environmental_interaction/S1_Table.xlsx⁠ from ⁠http://nematode.net⁠. S1_Table.xlsx was created by Dr. Bruce A Rosa for Martin et al. (2019). Permission to share this data was obtained from Dr. Bruce Rosa and Dr. Ivonne Martin.

References

Martin I, Uh H, Supali T, Mitreva M, Houwing-Duistermaat JJ (2019). “The mixed model for the analysis of a repeated-measurement multivariate count data.” Statistics in Medicine, 38(12), 2248–2268. doi:10.1002/sim.8101.

Scealy JL, Wood ATA (2023). “Score matching for compositional distributions.” Journal of the American Statistical Association, 118(543), 1811–1823. doi:10.1080/01621459.2021.2016422.

Wiria AE, Prasetyani MA, Hamid F, Wammes LJ, Lell B, Ariawan I, Uh HW, Wibowo H, Djuardi Y, Wahyuni S, Sutanto I, May L, Luty AJ, Verweij JJ, Sartono E, Yazdanbakhsh M, Supali T (2010). “Does treatment of intestinal helminth infections influence malaria?” BMC Infectious Diseases, 10, 77. doi:10.1186/1471-2334-10-77.

Description

Build new tapes (i.e ADFun objects) from existing tapes, including differentiation, swapping independent variables and dynamic parameters, and Jacobian determinants.

Usage

tapeJacobian(tape)

tapeHessian(tape)

tapeGradOffset(tape)

tapeLogJacDet(tape)

tapeSwap(tape)

Arguments

Details

The information in the fields xtape and dyntape of tape are used to perform the taping.

tapeJacobian

The returned vector is ordered with the range elements iterating fastest, then the domain elements. See https://cppad.readthedocs.io/latest/Jacobian.html.

tapeHessian

Suppose the function represented by tape maps from $d$-dimensional space to $1$-dimensional space, then the first $d$ elements of the vector is the gradient of the partial derivative with respect to the first dimension of the function's domain. The next $d$ elements of the vector is the gradient of the partial derivative of the second dimension of the function's domain. The Hessian as a matrix, can be obtained by using matrix() with ncol = d.

tapeGradOffset

A quadratic function can be written as

$f(x;\theta) = \frac{1}{2} x^T W(\theta) x + b(\theta)^Tx + c,$

where the vector $x$ is the independent variable of tape and the vector $\theta$ is the dynamic parameter vector of tape. The gradient of $f(x; \theta)$ with respect to $x$ is

$\Delta f(x; \theta) = \frac{1}{2}(W(\theta) + W(\theta)^T)x + b(\theta).$

The Hessian is

$H f(x; \theta) = \frac{1}{2}(W(\theta) + W(\theta)^T),$

which does not depend on $x$, so the gradient of the function can be rewritten as

$\Delta f(x;\theta) = H f(x; \theta) x + b(\theta)^T.$

The tape calculates $b(\theta)$ as

$b(\theta) = \Delta f(x;\theta) - H f(x; \theta) x,$

which does not depend on $x$.

Value

An ADFun object.

Functions

• tapeJacobian: Tape the Jacobian of a tape. The resulting tape returns the Jacobian as a vector.

• tapeHessian: Tape the Hessian of a tape. The resulting tape returns the Jacobian as a vector (see https://cppad.readthedocs.io/latest/Hessian.html).

• tapeGradOffset: A quadratic function can be written as

  $f(x;\theta) = \frac{1}{2} x^T W(\theta) x + b(\theta)^Tx + c.$

  The function tapeGradOffset creates a tape of $b(\theta)$ where $\theta$ is the independent variable.

• tapeLogJacDet: Creates a tape of the log of the Jacobian determinant of a function taped in tape. The dimensions of the domain (length of independent variable) and range (length of output variable) of tape must be equal for computation of the determinant.

• tapeSwap: Convert an ADFun so that the independent variables become dynamic parameters and the dynamic parameters become independent variables.

See Also

ADFun

Other tape builders: buildsmdtape()

Examples

tapes <- buildsmdtape("sph", "identity", "sph", "vMF",
              ytape = rep(1, 3)/sqrt(3),
              usertheta = rep(NA, 3)
              ) 
tapeJacobian(tapes$smdtape)
tapeHessian(tapes$smdtape)
tapeLogJacDet(tapeJacobian(tapes$smdtape))
tapeSwap(tapes$smdtape)

Description

Estimates the parameters of the Polynomially-Tilted Pairwise Interaction (PPI) model (Scealy and Wood 2023) for compositional data. By default ppi() uses cppad_closed() to find estimate. For many situations a hard-coded implementation of the score matching estimator is also available.

For a given parameter vector evalparam, ppi_smvalues() computes the score matching discrepancy, the gradient and the Hessian of the score matching discrepancy (see smvalues()) and the gradient offset of the score matching discrepancy (see quadratictape_parts() and tapeGradOffset()).

Usage

ppi(
  Y,
  paramvec = NULL,
  trans,
  method = "closed",
  w = rep(1, nrow(Y)),
  constrainbeta = FALSE,
  bdryw = "ones",
  acut = NULL,
  bdrythreshold = 1e-10,
  shiftsize = bdrythreshold,
  approxorder = 10,
  control = list(tol = 1e-15, checkgrad = TRUE),
  paramvec_start = NULL
)

ppi_smvalues(
  Y,
  paramvec = NULL,
  evalparam,
  trans,
  method = "closed",
  w = rep(1, nrow(Y)),
  bdryw = "ones",
  acut = NULL,
  bdrythreshold = 1e-10,
  shiftsize = bdrythreshold,
  approxorder = 10,
  average = TRUE
)

Arguments

Details

Estimation may be performed via transformation of the measure in Hyvärinen divergence from Euclidean space to the simplex (inverse of the additive log ratio transform), from a hyperplane to the simplex (inverse of the centred log ratio transform), from the positive quadrant of the sphere to the simplex (inverse of the square root transform), or without any transformation. In the latter two situations there is a boundary and weighted Hyvärinen divergence (Equation 7, Scealy and Wood 2023) is used. Properties of the estimator using the square root transform were studied by Scealy and Wood (2023). Properties of the estimator using the additive log ratio transform were studied by Scealy et al. (2024).

There are three boundary weight functions available:

• The function "ones" applies no weights and should be used whenever the manifold does not have a boundary.

• The function "minsq" is the minima-based boundary weight function for the PPI model (Equation 12, Scealy and Wood 2023)

  $\tilde{h}(z)^2 = \min(z_1^2, z_2^2, ..., z_p^2, a_c^2).$

  where $z$ is a point in the positive orthant of the p-dimensional unit sphere and $z_j$ is the jth component of z.

• The function "prodsq" is the product-based (Equation 9, Scealy and Wood 2023)

  $\tilde{h}(z)^2 = \min(\prod_{j=1}^{p} z_j^2, a_c^2).$

  where $z$ is a point in the positive orthant of the p-dimensional unit sphere and $z_j$ is the jth component of z.

Scealy and Wood (Theorem 1, Scealy and Wood 2023) prove that minimising the weighted Hyvärinen Divergence is equivalent to minimising $\psi(f, f_0)$ (See scorematchingtheory) when the boundary weight function is smooth or for the functions "minsq" and "prodsq" above when the manifold is the simplex or positive orthant of a sphere.

Hard-coded estimators are available for the following situations:

• Square root transformation ("sqrt") with the "minsq" boundary weight function:

  • full parameter vector (paramvec not provided)

  • paramvec fixes only the final element of $\beta$

  • paramvec fixes all elements of $\beta$

  • paramvec fixes $b_L = 0$ and provides fixed values of $\beta$

  • paramvec fixes $A_L=0$ and $b_L=0$, leaving $\beta$ to be fitted.

• Square root transformation ("sqrt") with the "prodsq" boundary weight function:

  • paramvec fixes all elements of $\beta$

  • paramvec fixes $b_L = 0$ and provides fixed values of $\beta$

  • paramvec fixes $A_L=0$ and $b_L=0$, leaving $\beta$ to be fitted.

• The additive log ratio transformation ("alr") using the final component on the denominator, with $b_L=0$ and fixed final component of $\beta$.

Value

ppi() returns: A list of est, SE and info.

• est contains the estimates in vector form, paramvec, and as $A_L$, $b_L$ and $\beta$.

• SE contains estimates of the standard errors if computed. See cppad_closed().

• info contains a variety of information about the model fitting procedure and results.

ppi_smvalues() returns a list of

• obj the score matching discrepancy value

• grad the gradient of the score matching discrepancy

• hess the Hessian of the score matching discrepancy

• offset gradient offset (see quadratictape_parts())

PPI Model

The PPI model density is proportional to

$\exp(z_L^TA_Lz_L + b_L^Tz_L)\prod_{i=1}^p z_i^{\beta_i},$

where $p$ is the dimension of a compositional measurement $z$, and $z_L$ is $z$ without the final ($p$th) component. $A_L$ is a $p-1 \times p-1$ symmetric matrix that controls the covariance between components. $b_L$ is a $p-1$ vector that controls the location within the simplex. The $i$th component $\beta_i$ of $\beta$ controls the concentration of density when $z_i$ is close to zero: when $\beta_i \geq 0$ there is no concentration and $\beta_i$ is hard to identify; when $\beta_i < 0$ then the probability density of the PPI model increases unboundedly as $z_i$ approaches zero, with the increasing occuring more rapidly and sharply the closer $\beta_i$ is to $-1$.

References

Scealy JL, Hingee KL, Kent JT, Wood ATA (2024). “Robust score matching for compositional data.” Statistics and Computing, 34, 93. doi:10.1007/s11222-024-10412-w.

Scealy JL, Wood ATA (2023). “Score matching for compositional distributions.” Journal of the American Statistical Association, 118(543), 1811–1823. doi:10.1080/01621459.2021.2016422.

See Also

Other PPI model tools: dppi(), ppi_param_tools, ppi_robust(), rppi()

Examples

model <- rppi_egmodel(100)
estalr <- ppi(model$sample,
              paramvec = ppi_paramvec(betap = -0.5, p = ncol(model$sample)),
              trans = "alr")
estsqrt <- ppi(model$sample,
              trans = "sqrt",
              bdryw = "minsq", acut = 0.1)

Description

These functions help to quickly generate a set of Windham exponents for use in ppi_robust() or Windham(). Rows and columns of $A_L$ and $b_L$ corresponding to components with strong concentrations of probability mass near zero have non-zero constant tuning exponent, and all other elements have a tuning constant of zero. All elements of $\beta$ have a tuning exponent of zero.

The function ppi_cW_auto() automatically detects concentrations near zero by fitting a PPI distribution with $A_L=0$ and $b_L=0$ (i.e. a Dirichlet distribution) with the centred log-ratio transformation of the manifold.

Usage

ppi_cW(cW, ...)

ppi_cW_auto(cW, Y)

Arguments

Details

The Windham robustifying method involves weighting observations by a function of the proposed model density (Windham 1995). Scealy et al. (2024) found that only some of the tuning constants should be non-zero: the tuning exponents corresponding to $\beta$ should be zero to avoid infinite weights;and to improve efficiency any rows or columns of $A_L$ corresponding to components without concentrations of probability mass (i.e. outliers can't exist) should have exponents of zero. Scealy et al. (2024) set the remaining tuning exponents to a constant.

Value

A vector of the same length as the parameter vector of the PPI model. Elements of $A_L$ will have a value of cW if both their row and column component has probability mass concentrated near zero. Similarly, elements of $b_L$ will have a value of cW if their row corresponds to a component that has a probability mass concentrated near zero. All other elements are zero.

References

Scealy JL, Hingee KL, Kent JT, Wood ATA (2024). “Robust score matching for compositional data.” Statistics and Computing, 34, 93. doi:10.1007/s11222-024-10412-w.

Windham MP (1995). “Robustifying Model Fitting.” Journal of the Royal Statistical Society. Series B (Methodological), 57(3), 599–609. 2346159, http://www.jstor.org/stable/2346159.

Examples

Y <- rppi_egmodel(100)$sample
ppi_cW_auto(0.01, Y)
ppi_cW(0.01, TRUE, TRUE, FALSE)

Description

Computes a marginal moment matching estimator (Section 6.2, Scealy and Wood 2023), which assumes $\beta$ is a known vector with the same value in each element, and $b_L = 0$. Only $A_L$ is estimated.

Usage

ppi_mmmm(Y, ni, beta0, w = rep(1, nrow(Y)))

Arguments

Details

$\beta=\beta_0$ is fixed and not estimated. $b_L$ is fixed at zero. See (Section 6.2 and A.8 of Scealy and Wood 2023). The boundary weight function in the score matching discrepancy is the unthresholded product weight function

$h(z)^2 = \min\left(\prod_{j=1}^{p} z_j^2, a_c^2\right).$

Value

A vector of estimates for $A_L$ entries (diagonal and off diagonal).

References

Scealy JL, Wood ATA (2023). “Score matching for compositional distributions.” Journal of the American Statistical Association, 118(543), 1811–1823. doi:10.1080/01621459.2021.2016422.

Description

The default parameterisation of the PPI model is a symmetric covariance-like matrix $A_L$, a location-like vector $b_L$ and a set of Dirichlet exponents $\beta$. For p components, $A_L$ has p-1 rows, $b_L$ is a vector with p-1 elements and $\beta$ is a vector with p elements. For score matching estimation this form of the parameters must be converted into a single parameter vector using ppi_paramvec(). ppi_paramvec() also includes easy methods to set parameters to NA for estimation with ppi() (in ppi() the NA-valued elements are estimated and all other elements are fixed). The reverse of ppi_paramvec() is ppi_parammats(). An alternative parametrisation of the PPI model uses a single p by p matrix $A^*$ instead of $A_L$ and $b_L$, and for identifiability $A^*$ is such that $1^T A^* 1 = 0$ where $1=(1,1,...,1)$ and $0=(0,0,..., 0)$ (Scealy and Wood 2023). Convert between parametrisations using ppi_toAstar() and ppi_fromAstar().

Usage

ppi_paramvec(
  p = NULL,
  AL = NULL,
  bL = NULL,
  Astar = NULL,
  beta = NULL,
  betaL = NULL,
  betap = NULL
)

ppi_parammats(paramvec)

ppi_toAstar(AL, bL)

ppi_fromAstar(Astar)

Arguments

Details

ppi_paramvec() returns a vector starting with the diagonal elements of $A_L$, then the off-diagonal elements extracted by upper.tri() (which extracts elements of $A_L$ along each row, left to right, then top to bottom), then $b_L$, then $\beta$.

The Astar parametrisation rewrites the PPI density as proportional to

$\exp(z^TA^*z)\prod_{i=1}^p z_i^{\beta_i},$

where $A^*$ (Astar) is a $p$ by $p$ matrix. Because $z$ lies in the simplex (in particular $\sum z_i = 1$), the density is the same regardless of the value of $1^T A^* 1$=sum(Astar), where $1$ is the vector of ones. Thus $A_L$ and $b_L$ specify $A^*$ up to an additive factor. In the conversion ppi_toAstar(), $A^*$ is returned such that $1^T A^* 1 = 0$. NULL values or NA elements are not allowed for ppi_toAstar() and ppi_fromAstar().

Value

ppi_paramvec(): a vector of length $p + (p-1)(2 + (p-1)/2)$.

ppi_parammats(): A named list of $A_L$, $b_L$, and $\beta$.

ppi_toAstar(): The matrix $A^*$.

ppi_fromAstar(): A list of the matrix $A_L$, the vector $b_L$ and a discarded constant.

PPI Model

The PPI model density is proportional to

$\exp(z_L^TA_Lz_L + b_L^Tz_L)\prod_{i=1}^p z_i^{\beta_i},$

where $p$ is the dimension of a compositional measurement $z$, and $z_L$ is $z$ without the final ($p$th) component. $A_L$ is a $p-1 \times p-1$ symmetric matrix that controls the covariance between components. $b_L$ is a $p-1$ vector that controls the location within the simplex. The $i$th component $\beta_i$ of $\beta$ controls the concentration of density when $z_i$ is close to zero: when $\beta_i \geq 0$ there is no concentration and $\beta_i$ is hard to identify; when $\beta_i < 0$ then the probability density of the PPI model increases unboundedly as $z_i$ approaches zero, with the increasing occuring more rapidly and sharply the closer $\beta_i$ is to $-1$.

See Also

Other PPI model tools: dppi(), ppi_robust(), ppi(), rppi()

Examples

ppi_paramvec(AL = "diag", bL = 0, betap = -0.5, p = 3)
vec <- ppi_paramvec(AL = rsymmetricmatrix(2), beta = c(-0.8, -0.7, 0))
ppi_parammats(vec)
Astar <- rWishart(1, 6, diag(3))[,,1]
ppi_fromAstar(Astar)
ppi_toAstar(ppi_fromAstar(Astar)$AL, ppi_fromAstar(Astar)$bL)

Description

ppi_robust() uses Windham() and ppi() to estimate a PPI distribution robustly. There are many arguments to the ppi() function and we highly recommend testing your arguments on ppi() first before running ppi_robust().

ppi_robust_alrgengamma() performs the Windham robustification algorithm exactly as described in Scealy et al. (2024) for score matching via log-ratio transform of the PPI model with $b_L = 0$. This function calls the more general Windham() and ppi().

Usage

ppi_robust(Y, cW, ...)

ppi_robust_alrgengamma(
  Y,
  cW,
  ...,
  fpcontrol = list(Method = "Simple", ConvergenceMetricThreshold = 1e-10)
)

Arguments

Details

ppi_robust_alrgengamma(): must fit a PPI model via additive-log ratio transform of the simplex with $b_L=0$ fixed and the final element of $\beta$ fixed. The default convergence metric and threshold are different to the default for ppi_robust() to match the implementation in (Scealy et al. 2024): convergence is measured by the change in the first element of $\beta$, and convergence is reached when the change is smaller than 1E-6. Override this behaviour by specifying the elements ConvergenceMetric and ConvergenceMetricThreshold in a list passed as fpcontrol. Windham() is called with alternative_populationinverse = TRUE.

Value

A list:

• est The estimated parameters in vector form (paramvec) and as AL, bL and beta.

• SE "Not calculated." Returned for consistency with other estimators.

• info Information returned in the optim slot of Windham(). Includes the final weights in finalweights.

References

Scealy JL, Hingee KL, Kent JT, Wood ATA (2024). “Robust score matching for compositional data.” Statistics and Computing, 34, 93. doi:10.1007/s11222-024-10412-w.

See Also

Other PPI model tools: dppi(), ppi_param_tools, ppi(), rppi()

Other Windham robustness functions: Windham(), vMF_robust()

Examples

set.seed(7)
model <- rppi_egmodel(100)
estsqrt <- ppi_robust(model$sample,
  cW = ppi_cW_auto(0.01, model$sample),
  paramvec_start = model$theta,
  trans = "sqrt", bdryw = "minsq", acut = 0.1)
set.seed(14)
model <- rppi_egmodel(100)
ppi_robust_alrgengamma(model$sample,
   cW = ppi_cW_auto(0.01, model$sample),
   paramvec = ppi_paramvec(betap = -0.5, p = ncol(model$sample)))

Description

When the score matching discrepancy function is quadratic then the gradient of the score matching discrepancy function can be written using the Hessian and an offset term. This can be useful for solving for the situation when the gradient is zero. The Hessian and offset term are computed using CppAD tapes. Taylor approximation can be used for locations at removed singularities (i.e. where intermediate values are unbounded). quadratictape_parts() will error if testquadratic(tape) returns FALSE.

Usage

quadratictape_parts(tape, tmat, tcentres = NA * tmat, approxorder = 10)

Arguments

Details

A quadratic function can be written

$f(x; t) = \frac{1}{2} x^T W(t) x + b(t)^T x + c,$

where $t$ is considered a vector that is constant with respect to the differentiation. The Hessian of the function is with respect to $x$ is

$H f(x; t) = \frac{1}{2}(W(t) + W(t)^T).$

The gradient of the function with respect to $x$ can then be written

$\Delta f(x;t) = H f(x; t) x + b(t)^T x,$

where the Hessian and offset $b(t)$ depend only on $t$.

The functions here evaluate the Hessian and offset $b(t)$ for many values of $t$. Tapes of the Hessian and gradient offset are created using tapeHessian() and tapeGradOffset() respectively. These tapes are then evaluated for every row of tmat. When the corresponding tcentres row is not NA, then approximate (but very accurate) results are calculated using Taylor approximation around the location given by the row of tcentres.

For score matching $x$ is the set of model parameters and the vector $t$ is a (multivariate) measurement.

Value

A list of

• offset Array of offsets $b(t)$, each row corresponding to a row in tmat

• Hessian Array of vectorised $H f(x; t)$ (see tapeHessian()), each row corresponding to a row in tmat. For each row, obtain the Hessian in matrix format by using matrix(ncol = length(tape$xtape)).

See Also

Other tape evaluators: evaltape(), smvalues(), testquadratic()

Examples

u <- rep(1/3, 3)
smdtape <- buildsmdtape("sim", "sqrt", "sph", "ppi",
              ytape = u,
              usertheta = ppi_paramvec(p = 3),
              bdryw = "minsq", acut = 0.01,
              verbose = FALSE
              )$smdtape
quadratictape_parts(smdtape, 
  tmat = rbind(u, c(1/4, 1/4, 1/2)))

Description

Given parameters of the PPI model, generates independent samples.

Usage

rppi(n, ..., paramvec = NULL, maxden = 4, maxmemorysize = 1e+05)

rppi_egmodel(n, maxden = 4)

Arguments

Details

We recommend running rppi() a number of times to ensure the choice of maxden is good. rppi() will error when maxden is too low.

The simulation uses a rejection-sampling algorithm with Dirichlet proposal (Appendix A.1.3 Scealy and Wood 2023). Initially n Dirichlet proposals are generated. After rejection there are fewer samples remaining, say $n^*$. The ratio $n^*/n$ is used to guess the number of new Dirichlet proposals to generate until n samples of the PPI model are reached.

Advanced use: The number of Dirichlet proposals created at a time is limited such that the matrices storing the Dirchlet proposals are always smaller than maxmemorysize bytes (give or take a few bytes for wrapping). Larger maxmemorysize leads to faster simulation so long as maxmemorysize bytes are reliably contiguously available in RAM.

Value

A matrix with n rows and p columns. Each row is an independent draw from the specified PPI distribution.

rppi_egmodel returns a list:

• sample A matrix of the simulated samples (n rows)

• p The number of components of the model

• theta The PPI parameter vector

• AL The $A_L$ parameter matrix

• bL The $b_L$ parameter vector

• beta The $\beta$ parameter vector

Functions

• rppi_egmodel: Simulates the 3-component PPI model from (Section 2.3, Scealy and Wood 2023) and returns both simulations and model parameters.

References

Scealy JL, Wood ATA (2023). “Score matching for compositional distributions.” Journal of the American Statistical Association, 118(543), 1811–1823. doi:10.1080/01621459.2021.2016422.

See Also

Other PPI model tools: dppi(), ppi_param_tools, ppi_robust(), ppi()

Examples

beta0=c(-0.8, -0.8, -0.5)
AL = diag(nrow = 2)
bL = c(2, 3)
samp <- rppi(100,beta=beta0,AL=AL,bL=bL)
rppi_egmodel(1000)

Description

A simple function for generating a symmetric matrix for use in examples. The diagonal, and upper-triangular elements of the matrix are simulated independently from a uniform distribution. The lower-triangle of the output matrix is copied from the upper-triangle. These matrices do not represent the full range of possible symmetric matrices.

Usage

rsymmetricmatrix(p, min = 0, max = 1)

Arguments

Value

A p x p symmetric matrix.

Examples

rsymmetricmatrix(5)

Description

This package includes score matching estimators for particular distributions and a general capacity to implement additional score matching estimators. Score matching is a popular estimation technique when normalising constants for the proposed model are difficult to calculate or compute. Score matching was first developed by Hyvärinen (2005) and was further developed for subsets of Euclidean space (Hyvärinen 2007; Yu et al. 2019; Yu et al. 2020; Liu et al. 2019), Riemannian manifolds (Mardia et al. 2016; Mardia 2018), and Riemannian manifolds with boundary (Scealy and Wood 2023). In this help entry we briefly describe score matching estimation.

Score Matching in General

In the most general form (Riemannian manifolds with boundary) score matching minimises the weighted Hyvärinen divergence (Equation 7, Scealy and Wood 2023)

$\phi(f, f_0) = \frac{1}{2} \int_M f_0(z)h(z)^2 \left\lVert P(z)\Big(\nabla_z \log(f) - \nabla_z \log(f_0)\Big)\right\rVert^2 dM(z),$

where

• $M$ is the manifold, isometrically embedded in Euclidean space, and $dM(z)$ is the unnormalised uniform measure on $M$.

• $P(z)$ is the matrix that projects points onto the tangent space of the manifold at $z$, which is closely related to to Riemannian metric of $M$.

• $f_0$ is the density of the data-generating process, defined with respect to $dM(z)$.

• $f$ is the density of a posited model, again defined with respect to $dM(z)$.

• $h(z)$ is a function, termed the boundary weight function, that is zero on the boundary of $M$ and smooth (Section 3.2, Scealy and Wood 2023) or potentially piecewise smooth.

• $\nabla_z$ is the Euclidean gradient operator that differentiates with respect to $z$.

• $\lVert \cdot \rVert$ is the Euclidean norm.

Note that, because $P(z)$ is the projection matrix, $\left\lVert P(z)\Big(\nabla_z \log(f) - \nabla_z \log(f_0)\Big)\right\rVert^2$ is the natural inner product of the gradient of the log ratio of $f$ and $f_0$.

When the density functions $f$ and $f_0$ are smooth and positive inside $M$, and the boundary weight function is smooth or of particular forms for specific manifolds (Section 3.2, Scealy and Wood 2023), then minimising the weighted Hyvärinen divergence $\phi(f, f_0)$ is equivalent to minimising the score matching discrepancy (Theorem 1, Scealy and Wood 2023)

$\psi(f, f_0) = \int f_0(z)\big(A(z) + B(z) + C(z)\big)dM(z),$

where

$A(z) = \frac{1}{2} h(z)^2 \left(\nabla_z \log(f)\right)^T P(z) \left(\nabla_z \log(f)\right),$

$B(z) = h(z)^2\Delta_z\log(f),$

$C(z) = \left(\nabla_z h(z)^2)\right)^T P(z) \left(\nabla_z \log(f)\right),$

and $\Delta$ is the Laplacian operator. We term

$A(z) + B(z) + C(z)$

the score matching discrepancy function.

We suspect that (Theorem 1, Scealy and Wood 2023) holds more generally for nearly all realistic continuous and piecewise-smooth boundary weight functions, although no proof exists to our knowledge.

When $n$ independent observations from $f_0$ are available, the integration in $\psi(f, f_0)$ can be approximated by an average over the observations,

$\psi(f, f_0) \approx \hat\psi(f, f_0) = \frac{1}{n} \sum_{i = 1}^n A(z_i) + B(z_i) + C(z_i).$

If we parameterise a family of models $f_\theta$ according to a vector of parameters $\theta$, then the score matching estimate is the $\theta$ that minimises $\hat\psi(f_\theta, f_0)$. In general, the score matching estimate must be found via numerical optimisation techniques, such as in the function cppad_search(). However, when the family of models is a canonical exponential family then often $\hat\psi(f_\theta, f_0)$ and the score matching discrepancy function is a quadratic function of $\theta$ (Mardia 2018) and the minimum has a closed-form solution found by cppad_closed().

Note that when $M$ has a few or more dimensions, the calculations of $A(z)$, $B(z)$ and $C(z)$ can become cumbersome. This package uses CppAD to automatically compute $A(z)$, $B(z)$ and $C(z)$, and the quadratic simplification if it exists.

Transformations

Hyvärinen divergence $\phi(f, f_0)$ is sensitive to transformations of the measure $dM(z)$, including transforming the manifold. That is, transforming the manifold $M$ changes the divergence between distributions and changes the minimum of $\hat\psi(f_\theta, f_0)$. The transformation changes measure $dM(z)$, the divergence and the estimator but does not transform the data.

For example, many different transformations of the simplex (i.e. compositional data) are possible (Appendix A.3, Scealy et al. 2024). Hyvärinen divergences that use the sphere, obtained from the simplex by a square root, have different behaviour to Hyvärinen divergence using Euclidean space obtained from the simplex using logarithms (Scealy et al. 2024). The estimator for the latter does not apply logarithms to the observations, in fact the estimator involves only polynomials of the observed compositions (Scealy et al. 2024).

The variety of estimator behaviour available through different transformations was a major motivator for this package as each transformation has different $A(z)$, $B(z)$ and $C(z)$, and without automatic differentiation, implementation of the score matching estimator in each case would require a huge programming effort.

References

Hyvärinen A (2005). “Estimation of Non-Normalized Statistical Models by Score Matching.” Journal of Machine Learning Research, 6(24), 695–709. https://jmlr.org/papers/v6/hyvarinen05a.html.

Hyvärinen A (2007). “Some extensions of score matching.” Computational Statistics & Data Analysis, 51(5), 2499–2512. doi:10.1016/j.csda.2006.09.003.

Liu S, Kanamori T, Williams DJ (2019). “Estimating Density Models with Truncation Boundaries using Score Matching.” doi:10.48550/arXiv.1910.03834.

Mardia K (2018). “A New Estimation Methodology for Standard Directional Distributions.” In 2018 21st International Conference on Information Fusion (FUSION), 724–729. doi:10.23919/ICIF.2018.8455640.

Mardia KV, Kent JT, Laha AK (2016). “Score matching estimators for directional distributions.” doi:10.48550/arXiv.1604.08470.

Scealy JL, Hingee KL, Kent JT, Wood ATA (2024). “Robust score matching for compositional data.” Statistics and Computing, 34, 93. doi:10.1007/s11222-024-10412-w.

Scealy JL, Wood ATA (2023). “Score matching for compositional distributions.” Journal of the American Statistical Association, 118(543), 1811–1823. doi:10.1080/01621459.2021.2016422.

Yu S, Drton M, Shojaie A (2019). “Generalized Score Matching for Non-Negative Data.” Journal of Machine Learning Research, 20(76), 1–70. https://jmlr.org/papers/v20/18-278.html.

Yu S, Drton M, Shojaie A (2020). “Generalized Score Matching for General Domains.” doi:10.48550/arXiv.2009.11428.

Description

Computes a range of relevant information for investigating score matching estimators.

Usage

smvalues(smdtape, xmat, pmat, xcentres = NA * xmat, approxorder = 10)

smvalues_wsum(
  tape,
  xmat,
  pmat,
  w = NULL,
  xcentres = NA * xmat,
  approxorder = 10
)

Arguments

Details

Computes the score matching discrepancy function from scorematchingtheory or weighted sum of the score matching discrepancy function. The gradient and Hessian are returned as arrays of row-vectors with each row corresponding to a row in xmat and pmat. Convert a Hessian row-vector to a matrix using matrix(ncol = length(smdtape$xtape)).

Value

A list of

• obj the score matching discrepancy values

• grad the gradient of the score matching discrepancy

• hess the Hessian of the score matching discrepancy

See Also

Other tape evaluators: evaltape(), quadratictape_parts(), testquadratic()

Examples

m <- rppi_egmodel(100)
smdtape <- buildsmdtape("sim", "sqrt", "sph", "ppi",
              ytape = rep(1/m$p, m$p),
              usertheta = ppi_paramvec(beta = m$beta),
              bdryw = "minsq", acut = 0.01)$smdtape
smvalues(smdtape, xmat = m$sample, pmat = m$theta[1:5])
smvalues_wsum(smdtape, m$sample, m$theta[1:5])$grad/nrow(m$sample)

Description

Uses the CppAD parameter property and derivatives (via tapeJacobian()) to test whether the tape is quadratic.

Usage

testquadratic(tape, xmat = NULL, dynparammat = NULL, verbose = FALSE)

Arguments

Details

Uses the xtape and dyntape values stored in tape to create new tapes. A tape of the Hessian is obtained by applying tapeJacobian() twice, and then uses a CppAD property to test whether the Hessian is constant. A function of quadratic form should have constant Hessian.

If xmat and dynparammat are non-NULL then testquadratic() also checks the Jacobian of the Hessian at xmat and dynparammat values. For quadratic form functions the Jacobian of the Hessian should be zero.

Value

TRUE or FALSE

See Also

Other tape evaluators: evaltape(), quadratictape_parts(), smvalues()

Examples

tapes <- buildsmdtape(
   "sim", "sqrt", "sph",
   ll = "ppi",
   ytape = c(0.2, 0.3, 0.5),
   usertheta = ppi_paramvec(p = 3), 
   bdryw = "minsq",
   acut = 0.1,
   verbose = FALSE)

 testquadratic(tapes$smdtape)

Description

In general the normalising constant in von Mises Fisher distributions is hard to compute, so Mardia et al. (2016) suggested a hybrid method that uses maximum likelihood to estimate the mean direction and score matching for the concentration. We can also estimate all parameters using score matching (smfull method), although this estimator is likely to be less efficient than the hybrid estimator. On the circle the hybrid estimators were often nearly as efficient as maximum likelihood estimators (Mardia et al. 2016). For maximum likelihood estimators of the von Mises Fisher distribution, which all use approximations of the normalising constant, consider movMF::movMF().

Usage

vMF(Y, paramvec = NULL, method = "Mardia", w = rep(1, nrow(Y)))

Arguments

Details

The full score matching estimator (method = "smfull") estimates $\kappa \mu$. The hybrid estimator (method = "Mardia") estimates $\kappa$ and $\mu$ separately. Both use cppad_closed() for score matching estimation.

Value

A list of est, SE and info.

• est contains the estimates in vector form, paramvec, and with user friendly names k and m.

• SE contains estimates of the standard errors if computed. See cppad_closed().

• info contains a variety of information about the model fitting procedure and results.

von Mises Fisher Model

The von Mises Fisher density is proportional to

$\exp(\kappa \mu^T z),$

where $z$ is on a unit sphere, $\kappa$ is termed the concentration, and $\mu$ is the mean direction unit vector. The effect of the $\mu$ and $\kappa$ can be decoupled in a sense (p169, Mardia and Jupp 2000), allowing for estimating $\mu$ and $\kappa$ separately.

References

Mardia KV, Jupp PE (2000). Directional Statistics, Probability and Statistics. Wiley, Great Britain. ISBN 0-471-95333-4.

Mardia KV, Kent JT, Laha AK (2016). “Score matching estimators for directional distributions.” doi:10.48550/arXiv.1604.08470.

See Also

Other directional model estimators: Bingham(), FB(), vMF_robust()

Examples

if (requireNamespace("movMF")){
  Y <- movMF::rmovMF(1000, 100 * c(1, 1) / sqrt(2))
  movMF::movMF(Y, 1) #maximum likelihood estimate
} else {
  Y <- matrix(rnorm(1000 * 2, sd = 0.01), ncol = 2)
  Y <- Y / sqrt(rowSums(Y^2))
}
vMF(Y, method = "smfull")
vMF(Y, method = "Mardia")
vMF(Y, method = "hybrid")

Description

Robust estimation for von Mises Fisher distribution using Windham().

Usage

vMF_robust(Y, cW, ...)

Arguments