empirical p-value function

.github/workflows/build-macos.yml

@@ -5,7 +5,7 @@ on: [push, workflow_dispatch]
 jobs:
   build:
 
-    runs-on: macos-11
+    runs-on: macos-15
 
     steps:
     - name: Checkout

.github/workflows/build-ubuntu.yml

@@ -5,7 +5,7 @@ on: [push]
 jobs:
   build:
 
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-24.04
 
     steps:
     - name: Checkout

CMakeLists.txt

@@ -1,5 +1,5 @@
 # package config
-cmake_minimum_required(VERSION 3.11.0)
+cmake_minimum_required(VERSION 3.24.0)
 project(
   LatAnalyze
   VERSION 3.6
@@ -41,18 +41,17 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/${CMAKE_INSTALL_LIBDIR}")
 
 # fetch and create symbolic link to Eigen headers
-set(FETCHCONTENT_BASE_DIR ${CMAKE_BINARY_DIR}/deps)
-FetchContent_Declare(
-  Eigen3
-  GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
-  GIT_TAG 3.4.0)
-FetchContent_GetProperties(Eigen3)
-if(NOT eigen3_POPULATED)
-  FetchContent_Populate(Eigen3)
-  message(STATUS "Eigen3 fetched")
+set(EIGEN_URL "https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz")
+set(EIGEN_TGZ "${CMAKE_BINARY_DIR}/eigen-3.4.0.tar.gz")
+set(EIGEN_DIR "${CMAKE_BINARY_DIR}/deps/eigen-src")
+if(NOT EXISTS "${EIGEN_DIR}/Eigen/Dense")
+    file(DOWNLOAD "${EIGEN_URL}" "${EIGEN_TGZ}" SHOW_PROGRESS)
+    file(ARCHIVE_EXTRACT INPUT "${EIGEN_TGZ}" DESTINATION "${CMAKE_BINARY_DIR}/deps")
+    file(GLOB EIGEN_EXTRACTED_DIR "${CMAKE_BINARY_DIR}/deps/eigen-*")
+    file(RENAME "${EIGEN_EXTRACTED_DIR}" "${EIGEN_DIR}")
+    message(STATUS "Eigen3 fetched")
 endif()
-file(CREATE_LINK ${eigen3_SOURCE_DIR}/Eigen
-     ${CMAKE_SOURCE_DIR}/lib/LatAnalyze/Eigen SYMBOLIC)
+file(CREATE_LINK ${EIGEN_DIR}/Eigen ${CMAKE_SOURCE_DIR}/lib/LatAnalyze/Eigen SYMBOLIC)
 
 # dependencies
 find_package(Threads REQUIRED)

ci-scripts/install-gsl.sh

@@ -1,6 +1,6 @@
 #!/usr/bin/env bash
 
-name='gsl-2.6'
+name='gsl-2.8'
 
 if (($# != 2)); then
   echo "usage: $(basename "$0") <prefix> <ntasks>" 1>&2
@@ -15,7 +15,7 @@ mkdir -p "${prefix}"
 cd "${prefix}"
 prefix=$(pwd -P)
 cd "${init_dir}/local/build"
-wget http://ftpmirror.gnu.org/gsl/${name}.tar.gz
+wget https://mirror.ibcp.fr/pub/gnu/gsl/${name}.tar.gz
 tar -xzvf ${name}.tar.gz
 mkdir -p ${name}/build
 cd ${name}/build

examples/exRand.cpp

@@ -32,9 +32,10 @@ int main(void)
     for (double s = 1.; s < 5.; ++s)
     {
         auto ci = h.confidenceInterval(s);
-        
+
         cout << static_cast<int>(s) << " sigma(s) interval= [";
         cout << ci.first << ", " << ci.second << "]" << endl;
+        cout << "P(X > " << s << ") = " << h.pValue(s) + 1. - h.pValue(-s) << endl;
     }
     p << PlotHistogram(h);
     p << PlotFunction(compile("return exp(-x_0^2/2)/sqrt(2*pi);", 1), -5., 5.);

lib/CMakeLists.txt

@@ -88,10 +88,12 @@ install(
   FILE LatAnalyzeTargets.cmake
   NAMESPACE LatAnalyze::
   DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/LatAnalyze)
+set(INCLUDE_INSTALL_DIR "${CMAKE_INSTALL_INCLUDEDIR}")
 configure_package_config_file(
   ${CMAKE_CURRENT_SOURCE_DIR}/Config.cmake.in
   "${CMAKE_CURRENT_BINARY_DIR}/LatAnalyzeConfig.cmake"
-  INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/LatAnalyze)
+  INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/LatAnalyze
+  PATH_VARS INCLUDE_INSTALL_DIR)
 write_basic_package_version_file(
   "${CMAKE_CURRENT_BINARY_DIR}/LatAnalyzeConfigVersion.cmake"
   VERSION "${PROJECT_VERSION}"

lib/Config.cmake.in

@@ -25,4 +25,6 @@ endforeach()
 
 include("${CMAKE_CURRENT_LIST_DIR}/LatAnalyzeTargets.cmake")
 
+set(@PROJECT_NAME@_INCLUDE_DIRS "@PACKAGE_INCLUDE_INSTALL_DIR@")
+
 check_required_components(LatAnalyze)

lib/LatAnalyze/Core/Math.cpp

@@ -166,5 +166,17 @@ auto chi2CcdfVecFunc = [](const double arg[2])
     return gsl_cdf_chisq_Q(arg[0], arg[1]);
 };
 
+auto hotellingT2PValueVecFunc = [](const double arg[3])
+{
+    double T2 = arg[0];
+    double n  = arg[1];
+    double p  = arg[2];
+    double F = (n - p) / (p * (n - 1)) * T2;
+    double p_value = 1.0 - gsl_cdf_fdist_P(F, p, n - p);
+
+    return p_value;
+};
+
 DoubleFunction MATH_NAMESPACE::chi2PValue(chi2PValueVecFunc, 2);
 DoubleFunction MATH_NAMESPACE::chi2Ccdf(chi2CcdfVecFunc, 2);
+DoubleFunction MATH_NAMESPACE::hotellingT2PValue(hotellingT2PValueVecFunc, 3);

lib/LatAnalyze/Core/Math.hpp

@@ -160,6 +160,7 @@ namespace MATH_NAMESPACE
 {
     extern DoubleFunction chi2PValue;
     extern DoubleFunction chi2Ccdf;
+    extern DoubleFunction hotellingT2PValue;
 }
 
 END_LATAN_NAMESPACE

lib/LatAnalyze/Core/Plot.cpp

@@ -511,6 +511,57 @@ PlotImpulses::PlotImpulses(const DVec &x, const DVec &y)
     setCommand("'" + tmpFileName + "' u 1:2 w impulses");
 }
 
+// PlotSteps constructor ////////////////////////////////////////////////////
+PlotSteps::PlotSteps(const DVec &x, const DVec &y)
+{
+    if (x.rows() != y.rows())
+    {
+        LATAN_ERROR(Size, "x and y vector does not have the same size");
+    }
+
+    DMat   d(x.rows(), 2);
+    string tmpFileName;
+
+    for (Index i = 0; i < x.rows(); ++i)
+    {
+        d(i, 0) = x(i);
+        d(i, 1) = y(i);
+    }
+    tmpFileName = dumpToTmpFile(d);
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' u 1:2 w steps");
+}
+
+// PlotGrid constructor ////////////////////////////////////////////////////////
+PlotGrid::PlotGrid(const DVec &x, const DVec &y, const DMat &value)
+{
+    if (x.size() != value.rows())
+    {        
+        LATAN_ERROR(Size, "x vector does not have the same size as value matrix rows");
+    }
+    if (y.size() != value.cols())
+    {        
+        LATAN_ERROR(Size, "y vector does not have the same size as value matrix columns");
+    } 
+    if (value.rows() < 2 || value.cols() < 2)
+    {
+        LATAN_ERROR(Size, "value matrix must have at least 2 rows and 2 columns");
+    }
+    DMat   d(value.cols()+1, value.rows()+1);
+    string tmpFileName;
+    d(0,0) = value.cols();
+    d.row(0).tail(value.cols()) = x;
+    d.col(0).tail(value.rows()) = y;
+    for (Index i = 0; i < value.rows(); ++i)
+    for (Index j = 0; j < value.cols(); ++j)
+    {
+        d(i+1, j+1) = value(j, i);
+    }
+    tmpFileName = dumpToTmpFile(d);
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' nonuniform matrix w image");
+}
+
 // PlotMatrixNoRange constructor ///////////////////////////////////////////////
 PlotMatrixNoRange::PlotMatrixNoRange(const DMat &m)
 {

lib/LatAnalyze/Core/Plot.hpp

@@ -191,6 +191,15 @@ public:
     virtual ~PlotHistogram(void) = default;
 };
 
+class PlotSteps: public PlotObject
+{
+public:
+    // constructor
+    PlotSteps(const DVec &x, const DVec &y);
+    // destructor
+    virtual ~PlotSteps(void) = default;
+};
+
 class PlotImpulses: public PlotObject
 {
 public:
@@ -200,6 +209,15 @@ public:
     virtual ~PlotImpulses(void) = default;
 };
 
+class PlotGrid: public PlotObject
+{
+public:
+    // constructor
+    PlotGrid(const DVec &x, const DVec &y, const DMat &value);
+    // destructor
+    virtual ~PlotGrid(void) = default;
+};  
+
 class PlotMatrixNoRange: public PlotObject
 {
 public:

lib/LatAnalyze/Numerical/DWT.cpp

@@ -38,6 +38,10 @@ void filterConvolution(MatType &out, const MatType &data,
 {
     Index n = data.rows(), nf = n*filter.size();
 
+    if (&out == &data)
+    {
+        LATAN_ERROR(Argument, "filter convolution does not support in-place operation");
+    }
     out.resizeLike(data);
     out.fill(0.);
     for (unsigned int i = 0; i < filter.size(); ++i)

lib/LatAnalyze/Statistics/Histogram.cpp

@@ -162,6 +162,20 @@ double Histogram::operator()(const double x) const
     return (*this)[static_cast<Index>(i)];
 }
 
+// p-value P(x > x0) ///////////////////////////////////////////////////////////
+double Histogram::pValue(const double x0) const
+{
+    Index n = data_.size();
+    double count = 0.;
+
+    FOR_STAT_ARRAY(data_, s)
+    {
+        count += (data_[s] > x0) ? 1. : 0.;
+    }
+    
+    return count/n;
+} 
+
 // percentiles & confidence interval ///////////////////////////////////////////
 double Histogram::percentile(const double p) const
 {

lib/LatAnalyze/Statistics/Histogram.hpp

@@ -54,6 +54,8 @@ public:
     double                    getX(const Index i) const;
     double                    operator[](const Index i) const;
     double                    operator()(const double x) const;
+    // p-value P(x > x0)
+    double                    pValue(const double x0) const;
     // percentiles & confidence interval
     double                    percentile(const double p) const;
     double                    median(void) const;