HistDAWass

(Histogram-valued Data analysis using Wasserstein

metric)

What’s new in version 1.0.9

Changes and improvements with respect to version 1.0.8.

HistDAWass 1.0.9 is a maintenance release. It fixes a memory bug reported by the Rcpp team, introduces a consistent treatment of empty (zero-weight) bins, fixes several numerical bugs and makes the C++ distance computations considerably faster, and revises the graphics for ggplot2 4.x. Results on data without empty bins are numerically identical to version 1.0.8.

1. Rcpp compatibility fix (Rcpp PR #1513)

The next Rcpp release checks that a vector assigned into a matrix row has the same length as the row. With that check the example of WH_adaptive.kmeans() failed with:

Error: cannot assign a vector of length 2 to a target of length 3

Cause. In c_WH_ADPT_KMEANS_TOTALSSQ() (src/Utility_fun.cpp) the vector of cluster weights was allocated with k elements (the number of clusters) but assigned into a row of a matrix with vars columns (the number of variables). When k < vars, released Rcpp read past the end of the vector, so the weights of the remaining variables came from uninitialised memory.

Fix. Each row of the weight matrix is now filled directly with vars copies of the cluster weight. The function is used by WH_adaptive.kmeans() and by the adaptive Kohonen maps.

Verified by running all the examples and R CMD check --as-cran against both the released Rcpp and the Rcpp development branch of PR #1513 (version 1.1.2.4, with the length checks enabled).

Two further out-of-bounds accesses in the C++ code were fixed:

• whichLowerEQThanX() could read one element beyond the end of the vector;
• c_DISTA_M2() accessed the second element of a list of length one when the data had a single variable.

2. Empty bins (zero-weight bins)

A histogram may have bins with zero weight, i.e. two consecutive equal values in the cdf p. In 1.0.8 they were handled in different and inconsistent ways (a weight of 0.001 counts in data2hist(), 1e-14 in the C++ registration, 1e-13 in checkEmptyBins(), nothing in the constructor).

New rule (everywhere): each empty bin gets a weight of 1e-8 and the cdf is renormalised so that it starts at 0 and ends at 1. Bins are detected with a threshold of 0.5e-8, so that applying the rule twice does not change an already fixed distribution.

The rule is applied:

• when a distributionH or TdistributionH object is created (distributionH(), new("distributionH", ...)); mean and standard deviation are recomputed on the fixed cdf;
• in data2hist() (all algorithms, including "ManualBreaks", that did not handle empty bins at all);
• in checkEmptyBins(), which now also has a method for MatH objects (useful for objects created with previous versions or modified by direct slot assignment);
• in the C++ functions that register distributions on a common grid of cdf values (REGISTER2() and PREPARE_A_VEC_MAT()), so that objects created with older versions are handled correctly too.

Why it matters — a real bug, not only a numerical nicety. In 1.0.8, when a distribution had repeated values of p, the C++ code that puts several distributions on a common grid misaligned the quantiles. Distances, clusterings and other statistics on such data were wrong: on the China_Seas data the L2 Wasserstein distances computed by WH_MAT_DIST() had relative errors up to 7% with respect to a brute-force numerical integration of the squared difference of the quantile functions. In 1.0.9 the error is about 1e-5 (due only to the 1e-8 weights and to the brute-force grid), also for objects created with older versions.

Datasets. The shipped datasets China_Seas (345 distributions), China_Month (296) and Age_Pyramids_2014 (42) contained empty bins and have been fixed with the same rule. Results on these datasets differ slightly from 1.0.8 (and are now correct).

3. Other bug fixes

4. Performance

The C++ functions computing squared L2 Wasserstein distances between quantile functions (c_MM_L2_SQ_WASS_D(), c_DISTA_M(), c_DISTA_M2(), c_Fast_D_Mat(), c_STEP_2_ADA_KMEANS(), c_ComputeFast_L2_SQ_WASS_DMAT(), c_ComputeFast_L2_SQ_WASS_DMAT_simple()) have been rewritten:

• no temporary vectors are allocated in the inner loops (Rcpp sugar expressions were replaced by plain loops on pointers to the data);
• the quantiles and the mean of the prototypes are read from the S4 objects once, outside the loops, instead of once for each individual;
• in c_WH_ADPT_KMEANS_TOTALSSQ() the distances between each histogram and the general prototype, which do not depend on the cluster, are computed once instead of k times.

The new kernels perform exactly the same floating point operations in the same order, so results are bit-for-bit identical to 1.0.8 (checked on 37 functions of the package, on data without empty bins).

Median times on a MatH of 300 histograms x 4 variables (R 4.3.3, Linux):

5. Graphics

WH.MultiplePCA() has a new argument graph (default TRUE) to switch off the graphs of FactoMineR::MFA().

6. Tests and package infrastructure

• Documentation revised and regenerated with roxygen2 8: descriptive titles, one usage per method, accurate return values, related functions documented together (e.g. ?WH.var.covar covers sums of squares, covariances and correlations; ?"Arith-distributionH" the arithmetic operators), and an overview of the package in ?HistDAWass.

• Added a testthat test suite (tests/testthat) covering empty bins, the shipped datasets, compQ()/compP(), the validity checks, the adaptive k-means with k lower than the number of variables, and the agreement of WH_MAT_DIST() with pairwise WassSqDistH() and with a brute-force computation of the distance. The tests fail on 1.0.8 and pass on 1.0.9.

• DESCRIPTION: version 1.0.9, Suggests: testthat, URL and BugReports pointing to https://github.com/Airpino/HistDAWass; removed the fields generated by CRAN (Packaged, Repository, Date/Publication, MD5).

• Added NEWS.md, cran-comments.md, an RStudio project file and a .Rbuildignore.

7. Compatibility notes

• Results are unchanged for data without empty bins.
• Results change slightly (and become correct) for data with empty bins, including the China_Seas, China_Month and Age_Pyramids_2014 datasets.
• data2hist() with algo = "base" or "histogram" gives empty bins a weight of 1e-8 instead of the equivalent of 0.001 observations: differences in the cdf are of the order of 1e-6.
• Creating a distributionH with a decreasing cdf is now an error.
• Plot methods now draw the plot when called (also inside functions) and return the ggplot invisibly; WH.1d.PCA() draws on the current device.

8. Check results

R CMD check --as-cran: 0 errors, 0 warnings, only NOTEs related to the local build environment, with both the released Rcpp and the Rcpp development branch of PR #1513.

In this document we describe the main features of the HistDAWass package. The name is the acronym for Histogram-valued Data analysis using Wasserstein metric. The implemented classes and functions are related to the analysis of data tables containing histograms in each cell instead of the classical numeric values.

In this document we describe the main features of the HistDAWass package. The name is the acronym for Histogram-valued Data analysis using Wasserstein metric. The implemented classes and functions are related to the anlysis of data tables containing histograms in each cell instead of the classical numeric values.

What is the L2 Wasserstein metric?

given two probability density functions f and g, each one has a cumulative distribution function F and G and thei respectively quantile functions (the inverse of a cumulative distribution function) Q_f and Q_g. The L2 Wasserstein distance is

The implemented classes are those described in the following table

library(HistDAWass)
mydist=distributionH(x=c(0,1,2),p=c(0,0.3,1))

From raw data to histograms

data2hist functions

Basic statistics for a distributionH (A histogram)

• mean

  • the mean of a histogram

• standard deviation

  • the standard deviation of a histogram

• skewness

  • the third standardized moment of a histogram

• kurthosis

  • the fourth standardized momemt of a histogram

Basic statistics for a MatH (A matrix of histogrm-valued data)

• The average hisogram of a column
  • It is an average histogram that minimizes the sum of squared Wasserstein distances.
• The standard deviation of a variable
  • It is a number that measures the dispersion of a set of histograms.
• The covarince matrix of a MatH
  • It is a matrix that measures the covariances into a set of hitogram variables.
• The correlation matrix of a MatH
  • It is a matrix that measures the correlation into a set of hitogram variables.

Visualization > plot of a distributionH

plot of a MatH

plot of a HTS

Data Analysis methods

Clustering

• Kmeans

• Adaptive distance based Kmeans

• Fuzzy cmeans

• Fuzzy cmeans based on adaptive Wasserstein distances

• Kohonen batch self organizing maps

• Kohonen batch self organizing maps with Wasserstein adaptive distances

• Hierarchical clustering

Dimension reduction techniques

• Principal components analysis of a single histogram variable

• Principal components analysis of a set of histogram variables (using Multiple Factor Analysis)

Methods for Histogram time series

Smoothing

• Moving averages

• Exponential smoothing

Forecasting

• KNN prediction of histogram time series

Linear regression

A two component model for a linear regression using Least Square method