Skip to content

Commit

Permalink
Merge pull request #4 from x0range/development
Browse files Browse the repository at this point in the history
Fixed README and DESCRIPTION and documentation
  • Loading branch information
x0range authored Jul 23, 2020
2 parents 16e7138 + 1e10dd1 commit 3275b43
Show file tree
Hide file tree
Showing 6 changed files with 124 additions and 37 deletions.
25 changes: 15 additions & 10 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,19 +1,24 @@
Package: finity
Type: Package
Title: Test for finiteness of moments in a distribution
Version: 0.1
Date: 2019-11-14
Title: Test for Finiteness of Moments in a Distribution
Version: 0.1.4
Date: 2020-04-17
Authors@R: c(person("Torsten", "Heinrich", role = c("aut", "cre"),
email = "torsten.heinrich@oxfordmartin.ox.ac.uk"),
email = "torsten.heinrich@posteo.net"),
person("Julian", "Winkler", role = c("aut"),
email = "[email protected]"))
Author: Torsten Heinrich <[email protected]>,
Julian Winkler <[email protected]>
Maintainer: Torsten Heinrich <[email protected]>
Description: Tests whether a given moment of the distribution of a given
sample is finite or not.
Author: Torsten Heinrich [aut, cre], Julian Winkler [aut]
Maintainer: Torsten Heinrich <[email protected]>
Description: The purpose of this package is to tests whether a given
moment of the distribution of a given sample is finite or not. For
heavy-tailed distributions with tail exponent b, only moments of
order smaller than b are finite. Tail exponent and heavy-
tailedness are notoriously difficult to ascertain. But the
finiteness of moments (including fractional moments) can be tested
directly. This package does that following the test suggested by
Trapani (2016) <doi:10.1016/j.jeconom.2015.08.006>.
License: GPL (>= 2)
Imports: Rcpp (>= 1.0.3)
Imports: Rcpp (>= 1.0.3), stabledist (>= 0.7)
LinkingTo: Rcpp, RcppArmadillo, BH
RoxygenNote: 7.0.0
Encoding: UTF-8
11 changes: 11 additions & 0 deletions R/example.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
.example_finite_moment_test <- function() {

# Generate sample
rvs <- stabledist::rstable(100000, 1.9, 0.5, 1, 0, pm = 0)

# Perform test
result <- finite_moment_test(rvs, 2)

# Print results
message(paste("Test statistic:", result[1], "p-value:", result[2], "\n\n"))
}
58 changes: 58 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
@@ -1 +1,59 @@
# finity

Package for testing whether a given moment of the distribution of a given sample is finite or not.


# Details

For heavy-tailed distributions with tail exponent a, only moments of order < a are finite. The tail
index and heavy- tailedness are is notoriously difficult to ascertain. But the finiteness of moments
(including fractional moments) can be tested directly. This package does that following the test
suggested by Trapani (2016).


# Installing from CRAN

```install.packages(finity)```


# Usage: How to perform the finite moment test

```
library(stabledist)
library(finity)
# Generate test data
rvs <- rstable(10000000, 1.9, 0.5, 1, 0, pm = 0)
# Perform the test
result <- finite_moment_test(rvs, 2)
```


# Building and installing the package from github

```
Rscript -e 'setwd("finity");Rcpp::compileAttributes()'
R CMD build finity
R CMD INSTALL finity_0.1.1.tar.gz
```
In the install command, ```finity_0.1.1.tar.gz``` should be replaced with the correct filename of the archive (in particular with the correct version number).

```Rcpp::compileAttributes()``` can also be run from the R interpreter. The shell commands can alternatively be run from the R interpreter as ```system("R CMD build finity")``` etc.


# Building the pdf manual from the github package

```
Rscript -e 'setwd("finity");roxygen2::roxygenize();roxygen2::roxygenize()'
R CMD Rd2pdf finity
```

Alternatively, ```roxygen2::roxygenize()``` can be executed twice from the R interpreter before building the pdf with ```R CMD Rd2pdf finity```.

It is important to execute ```roxygenize()``` twice since the help files for the individual functions are only created the second time.


# Building as an R/CRAN package

Remove the ```examples``` directory and follow the steps for building (and installing) and building the pdf manual.
3 changes: 3 additions & 0 deletions examples/example_1.R
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,10 @@ sourceCpp("src/finite_moment_test.cpp")
# Performance test 1: Some heavy-tailed distribution
message("\nPerformance test 1: Large sample\n\n")

# Create test data
rvs <- rstable(50000000, 1.9, 0.5, 1, 0, pm = 0)

# Perform and time finite moment test
system.time(
{
result <- finite_moment_test(rvs, 2)
Expand Down
23 changes: 11 additions & 12 deletions man/finity-package.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -9,18 +9,18 @@
\packageDescription{finity}
}
\details{
The main function of the package, which performs the test for the finiteness of moments, is

The DESCRIPTION file:
\packageDESCRIPTION{finity}
\packageIndices{finity}
For heavy-tailed distributions with tail exponent \eqn{a}, only
moments of order \eqn{<a} are finite. The tail index and heavy-
tailedness are is notoriously difficult to ascertain. But
the finiteness of moments (including fractional moments) can
be tested directly. This package does that following the
test suggested by Trapani (2016).
finite_moment_test().

The main function is finite_moment_test().}
The test assumes a sample drawn from an unknown distrubution \eqn{F}. Given this sample, the finiteness or infiniteness of any moment of order \eqn{k} of distribution \eqn{F} can be ascertained. For this, the test follows a randomised testing procedure with artificial randomness. The absolute sample moment \eqn{\mu_k} of the desired order \eqn{k} of the sample is transformed into a test statistic, which follows a \eqn{\chi^2} distribution with one degree of freedom exactly if the moment of the same order \eqn{k} of \eqn{F} is not finite. The null hypothesis in the test is that the moment is infinite; the alternative is that it is finite.

It should be noted that while the moment of order \eqn{k} of \eqn{F} may be infinite, the sample moment \eqn{\mu_k} is always finite because the sample is of finite size. The sample moment will, however, diverge with growing sample size if the moment of the same order \eqn{k} of the original distribution \eqn{F} is not finite.

The test works as follows: A standard normal distribution is rescaled with \eqn{\sqrt{\exp(\mu_k)}}, yielding a normal distribution with mean 0 and either finite or infinite variance, depending on whether the hypothesis holds. For every observation of the resulting distribution, it is then tested, if the observation is located within an interval \eqn{[-u, u]}. The resulting binary quantity \eqn{\zeta} (0, or 1, true or false) follows a Bernoulli distribution with mean 1/2 exactly if the \eqn{k}th moment of \eqn{F} is infinite. Otherwise the mean is not 1/2. Sampling a number of different intervals characterized by different bounds \eqn{u} drawn from a distribution with finite support, the test aggregates over quantities \eqn{\zeta} such that the resulting test statistic follows a \eqn{\chi^2} exactly if \eqn{E(\zeta)=1/2}, i.e., if the \eqn{k}th moment of \eqn{F} is infinite.

Trapani (2016) offers some insights into the performance of the test and the impact its parameters have. These parameters are optional arguments of the testing function finite_moment_test().
}
\author{
\packageAuthor{finity}

Expand All @@ -38,9 +38,8 @@ https://github.com/x0range/finity
%~~ \code{\link[<pkg>:<pkg>-package]{<pkg>}} ~~
}
\examples{
library(stabledist)
# Generate sample
rvs <- rstable(50000000, 1.9, 0.5, 1, 0, pm = 0)
rvs <- stabledist::rstable(100000, 1.9, 0.5, 1, 0, pm = 0)
# Perform test
result <- finite_moment_test(rvs, 2)
# Print results
Expand Down
41 changes: 26 additions & 15 deletions src/finite_moment_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
// [[Rcpp::depends(BH)]]
// [[Rcpp::plugins("cpp11")]]

#define BOOST_DISABLE_ASSERTS

/* Following two lines to suppress deprecation warnings about
* integer_log2.hpp.
*/
Expand Down Expand Up @@ -135,8 +137,9 @@ arma::uvec build_zeta_unsorted(arma::vec xim, arma::vec u_series) { // This
return(zeta);
}

//' Returns Chi^2(1) percentile for test.
//' Chi^2(1) Percentile
//'
//' @description Returns the Chi^2(1) percentile for the test statistic.
//' @param value Chi^2(1) value (type: double).
//' @return Chi^2(1) percentile (type: double).
//' @examples
Expand All @@ -158,14 +161,14 @@ double get_chisq1_percentile(double value) {
}
}

//' Computes absolute moment of order k in sample of observations obs.
//'
//' Absolute Moment of Order k
//'
//' @description Computes the absolute moment of order k of a sample of observations.
//' @param obs Observations (type: armadillo numeric vector).
//' @param k Moment order (type: double)
//' @return Moment value (type: double)
//' @examples
//' library(stabledist)
//' rvs <- rstable(50000000, 1.9, 0.5, 1, 0, pm = 0)
//' rvs <- stabledist::rstable(100000, 1.9, 0.5, 1, 0, pm = 0)
//' absolute_moment <- compute_absolute_moment(rvs, 2)
//' @export
// [[Rcpp::export]]
Expand All @@ -191,21 +194,22 @@ double compute_absolute_moment(arma::vec obs, double k) {
return(mu);
}

//' Computes Trapani's (2016) finite moment test
//'
//' Finite Moment Test
//'
//' @description Computes Trapani's (2016) finite moment test for moment of order k of the distribution of a given the sample of observations obs. Knowledge of the identity of the distribution is not required. The null hypothesis is that the moment is infinite; the alternative is that it is finite. The function takes parameters of the test as optional arguments; some insights into the impact the choice of parameter values has are given in Trapani (2016).
//' @param obs Observations (type: armadillo numeric vector).
//' @param k Moment order (type: double)
//' @param r Artificial sample size (type: int). Default is N^0.8.
//' @param psi Pescaling moment (type: double). Must be <k. Default is 2.0.
//' @param u Sampling range width for sampling range [-u, u] (type: double) Default is 1.0.
//' @param force_random_variate_sample If True, draw random variates for xi and u_series. If False, use quantile function values from a regular percentile space grid. This represents the density function better. Defaiult is False.
//' @param ignore_errors Ignore errors caused by Inf and NaN results for too large absolute moments. If True, it will return test statistic=NA, pvalue=1. If False, it will stop with an error. Default is False. But normally this will indicate an infinite moment.
//' @param verbose If True, print detailed output for debugging. Default is False.
//' @param random_salting Salt number to be added to the random seed (type: int). This prevents identical random variate series if multiple instances are started and run in parallel. Default is 0.
//' @return Trapani's Theta test statistic (type: double).
//' @return Corresponding p-value (Chi^2(1) percentile) (type: double.
//' @return Corresponding p-value (Chi^2(1) percentile) (type: double).
//' @examples
//' library(stabledist)
//' rvs <- rstable(50000000, 1.9, 0.5, 1, 0, pm = 0)
//' rvs <- stabledist::rstable(100000, 1.9, 0.5, 1, 0, pm = 0)
//' result <- finite_moment_test(rvs, 2)
//' @export
// [[Rcpp::export]]
Expand All @@ -229,6 +233,9 @@ arma::vec finite_moment_test(arma::vec obs,
* force_random_variate_sample (bool): Draw random variates for xi and u_series. Default is using quantile
* function values from a regular percentile space grid. This represents
* the density function better.
* ignore_errors (bool): Ignore errors caused by Inf and NaN results for too large absolute moments. If True,
* it will return test statistic NA, pvalue 0. If False, it will stop with an error.
* Default is False. But normally this will indicate an infinite moment.
* verbose (bool): Print detailed output for debugging.
* random_salting (int): salt number to be added to the random seed. Prevents identical random variate series
* if multiple instances are started and run in parallel.
Expand Down Expand Up @@ -273,7 +280,7 @@ arma::vec finite_moment_test(arma::vec obs,
long double mu_psi = compute_absolute_moment(obs, psi);
// rescaling to mu* (Trapani 2016 eq 16)
if (verbose) Rcpp::Rcout << "mu is: " << mu << std::endl;
mu = mu / pow(mu_psi, k / psi);
mu = mu / pow(((long double)mu_psi), ((long double)(k / psi)));
if (verbose) Rcpp::Rcout << " ...rescaled to: " << mu << std::endl;
long double rescaling_factor_2 = pow(overloaded_std_norm_moment(psi), k / psi) / overloaded_std_norm_moment(k);
// rescaling to mu** (Trapani 2016 eq 17)
Expand All @@ -284,11 +291,13 @@ arma::vec finite_moment_test(arma::vec obs,
double exp_mu_half = long_exp_mu_half;
if (boost::math::isinf(exp_mu_half)) {
//Stop
Rcpp::Rcout << "Error: Absolute moment is too large. exp(mu/2) cannot be represented as double, which we need to do in the armadillo vector." << std::endl;
Rcpp::Rcout << " However, at this kind of value, you can safely assume that your moment is infinite. Any Trapani test would return p=0 for H0 (moment finite)." << std::endl;
Rcpp::Rcout << "Error: Rescaled absolute moment is too large. exp(mu/2) cannot be represented as double, which we need to do in the armadillo vector." << std::endl;
Rcpp::Rcout << " However, at this kind of value, you can safely assume that your moment is infinite. The test would return p=0 for H0 (moment finite)." << std::endl;
Rcpp::Rcout << " Absolute moment mu was in long double: " << long_exp_mu_half << ". In double it was: " << exp_mu_half << std::endl;
if (ignore_errors) {
arma::vec return_values = {NAN, 1.0};
//arma::vec return_values = {NAN, 1.0};
arma::vec return_values(2);
return_values << NAN << 1.0;
return(return_values);
} else {
Rcpp::stop("Infinite value detected!\n");
Expand Down Expand Up @@ -398,7 +407,9 @@ arma::vec finite_moment_test(arma::vec obs,
}

// Return
arma::vec return_values = {Theta, chisq1_percentile};
//arma::vec return_values = {Theta, chisq1_percentile};
arma::vec return_values(2);
return_values << Theta << chisq1_percentile;
return(return_values);
}

0 comments on commit 3275b43

Please sign in to comment.