kthohr / stats Goto Github PK

I cannot get this to compile. The error is stats_incl/rand/rnorm.ipp:114:27: error: no matching function for call to 'rnorm'. What am I doing wrong?

Thanks!

#define STATS_ENABLE_BLAZE_WRAPPERS
#include <stats.hpp>

int main(int argc, const char* argv[])
{
    stats::rand_engine_t engine;
    unsigned P = 2;

    blaze::DynamicMatrix<double> mu = blaze::ZeroMatrix<double>(P, 1);
    blaze::DynamicMatrix<double> Sigma = blaze::IdentityMatrix<double>(P);

    // All of these make the same compile error.  Which, if any, is the correct idiom?
    stats::rmvnorm(mu, Sigma, engine);
    stats::rmvnorm<blaze::DynamicMatrix<double>>(mu, Sigma, engine);
    stats::rmvnorm(P, mu, Sigma, engine);
}

Hi,
Thanks for the great package.
I tried to compute p-value for chi-squared distribution using the following code:
stats::pchisq(degrees_freedom,x_val)
Comparing to chi-squared table, the values are usually off more than 0.01.
Could you help out?
Thanks,
Bruno

Calling stats:pf(Inf, 2,2) or stats:df(Inf, 2,2) would return NaNs instead of the expected 1.0 and 0.0, respectively.

It'd be great to have some tags for the release versions like 2.2.0.

PR #28 fixes this issue.

Hi there!

Thanks for the stats codes, very helpful. I've found that the uniform distributions density function (dunif) returns a non-zero value even when the random variable is outside of the bounds. For example, consider the uniform distribution unif(0.1, 0.2). A random variable with value 0.15 will have density 10, and if the RV value is 0.201 its density should be 0, theoretically. I found that the current implementation of dunif might not meet this criterion.

Best,

Hi Keith,
Thanks for sharing your library.

Changing the define at the bottom of statslib_options.hpp,
from:
#define __stats_pointer_settings__ __restrict__

to:
#define __stats_pointer_settings__ __restrict

will make the library compatible with MSVC as well as GCC/Clang.

Thanks again,
Rafael.

It is not obvious to me what the syntax is for passing a vector of probabilities to stats::qnorm(). How would I replicate this R code.

a = qnorm(p= c(0.01, 0.20, 0.50, 0.80, 0.99))

I was able to pass single values to stats::qnorm, but got errors when I would attempt to pass a vector.

I am not quite sure how to use the multivariate normal distribution function using Eigen. The signature of the function appears to expect X, mu_par, and Sigma_par to all have the same type; however, in Eigen there is no generic "matrix" class (as there appears to be for Armadillo, according to the documentation) but each value has a specified size, i.e.

Eigen::Matrix<double, 3, 1> X(1, 1, 1);
Eigen::Matrix<double, 3, 1> mean(0, 0, 0);
Eigen::Matrix<double, 3, 3> Sigma = Eigen::Matrix3d::Identity();

// Or, more succinctly
// Eigen::Vector3d X(1, 1, 1);
// Eigen::Vector3d mean(0, 0, 0);
// Eigen::Matrix3d Sigma = Eigen::Matrix3d::Identity();

p = stats::dmvnorm(X, mean, Sigma);

This fails to compile with the following:

main.cpp:11:16: error: no matching function for call to 'dmvnorm'
  std::cout << stats::dmvnorm(X, mean, Sigma) << std::endl;
               ^~~~~~~~~~~~~~
./include/dens/dmvnorm.ipp:28:1: note: candidate template ignored: deduced conflicting types for parameter 'mT' ('Matrix<[2 * ...], 1, [2 * ...], 1>' vs. 'Matrix<[2 * ...], 3, [2 * ...], 3>')
dmvnorm(const mT& X, const mT& mu_par, const mT& Sigma_par, bool log_form)
^
1 error generated.

How do I go about using the multivariate functions with Eigen?

EDIT: I am able to get the above example to compile using the Eigen::MatrixXd classes which does serve as a "generic" matrix class. This will work for me. However, out of curiosity is there a way to use the stats library with compile-time Eigen sizes (like in my example)? Eigen recommends using fixed sizes known at compile time for small matrices and I would like to take advantage of those optimizations if possible.

The following program produces a segmentation fault for me:

#include "include/stats.hpp"
#include <iostream>

main() {
        std::cout << stats::qbeta(0.0022, 45.5, 80.5) << std::endl;
}

(Compile with e.g. g++ example.cc)

Valgrind reports that the fault is caused by a stack overflow in gcem::sqrt_recur. As a comparison, R gives the following:

> qbeta(0.0022, 45.5, 80.5)
[1] 0.2461071

This code here should return 1.0, it returns Inf
std::cout << stats::ppois(100, 5.0) << "\n";

Is there any support for Parameter estimation for GEV distribution?

Like gevfit() in MATLAB.

Thanks for sharing your work!

It would be interesting to document the library and provide examples. At least list supported distribution functions and features.

Thanks again!

DJuego

Hi! First of all, I'd like to thank for the great library! I'm playing around with the library and got following issue. My sample code is:

#include <iostream>
#include "stats.hpp"
#include <Eigen/Dense>

#define STATS_USE_EIGEN

int main() {

    Eigen::MatrixXd gamma_rvs = stats::rgamma<Eigen::MatrixXd>(100,50,3.0,2.0);
    std::cout << gamma_rvs << std::endl;

    return 0;
}

And compilation gave me:

/usr/local/bin/cmake --build /Users/alex/Development/stats_tests/cmake-build-debug --target stats_tests -- -j 4
Scanning dependencies of target stats_tests
[ 50%] Building CXX object CMakeFiles/stats_tests.dir/main.cpp.o
In file included from /Users/alex/Development/stats_tests/include/rand/rand.hpp:27:0,
from /Users/alex/Development/stats_tests/include/stats.hpp:37,
from /Users/alex/Development/stats_tests/main.cpp:2:
/Users/alex/Development/stats_tests/include/rand/rgamma.hpp:42:4: warning: inline function 'mT stats::rgamma(stats::uint_t, stats::uint_t, eT, eT) [with mT = Eigen::Matrix<double, -1, -1>; eT = double; stats::uint_t = unsigned int]' used but never defined
mT rgamma(const uint_t n, const uint_t k, const eT shape_par, const eT scale_par);
^~~~~~
[100%] Linking CXX executable stats_tests
Undefined symbols for architecture x86_64:
"Eigen::Matrix<double, -1, -1, 0, -1, -1> stats::rgamma<Eigen::Matrix<double, -1, -1, 0, -1, -1>, double>(unsigned int, unsigned int, double, double)", referenced from:
_main in main.cpp.o
ld: symbol(s) not found for architecture x86_64
collect2: error: ld returned 1 exit status
make[3]: *** [stats_tests] Error 1
make[2]: *** [CMakeFiles/stats_tests.dir/all] Error 2
make[1]: *** [CMakeFiles/stats_tests.dir/rule] Error 2
make: *** [stats_tests] Error 2

Could you please help me to figure out what's going on.
Thanks!
Alex

Hello,

Thanks for your work on this project. I'm trying to use the qpois function with different rates. The output is correct when rate=10 in the test cases, but if say the rate is 1135, the function returns all 66. See the test output below:

*** qpois: begin tests. ***

qpois(0.00): 66.00000. Success = 0
qpois(0.00): 66.00000. Success = 0
qpois(0.00): 66.00000. Success = 0
qpois(0.01): 66.00000. Success = 0
qpois(0.07): 66.00000. Success = 0
qpois(0.33): 66.00000. Success = 0
qpois(0.58): 66.00000. Success = 0
qpois(0.92): 66.00000. Success = 0
qpois(1.00): 66.00000. Success = 0

Thank you for any guidance on correcting this. Let me know if you need help debugging.

I'm working on a project where I need to keep backwards compatibility with GSL's behaviour in an application but use code under an Apache-2.0 license. I need to use the mt19937 with a particular seed rather than the mt19937_64 rng. Currently I believe the only way of doing this in the stats library is editing the stats_options.hpp file to use mt19937, set my own seed, and then ensuring I pass in a smaller type (i.e. floats rather than double) to functions I call to ensure the random number generator is used the same way, and generates the same sequence of numbers for testing. This works for most functions (although some, deep down, use larger types than I need and thus generates a difference sequence - but this isn't common). This does mean forcing downstream users to fetch and edit this options.hpp file, which will likely cause confusion.

If you currently have a way of overriding the random number generator without editing the stats_options.hpp file that'd be great to understand. If not, please let me know, and I'll add in a new option define, test it, and contribute it back as a PR. Please let me know how you'd like me to proceed.

It's a great library. I hope to use the OpenMP functionality soon too. Thanks for the hard work!

Does it contain Generalized extreme value distribution also?

E.g.

dataN = [5.44030563610939, 10.7071212794335, 6.98863965149785, 5.75855288152361, 4.35389980467031, 6.65717299986184,
         6.23022372943957, 4.61731971899753, 8.78179518233592, 5.68478054008583, 5.20860732492683, 4.21808687299822,
         8.64324644243904, 5.15451105578527, 7.03633883643414, 9.30334793245299, 5.40342918996923, 5.89134257860982,
         6.89237263218421, 5.00530206799751]

Matlab :
param = gevfit(dataN);

Hello. I'm trying to use Eigen wrapper. Every environment is okay except MSVC.

process : intel, 64bit
c++ compiler : clang++(15.0.7)
Generator : Visual Studio 17 2022
Eigen : 3.4.0

This is my code.

#define STATS_ENABLE_STDVEC_WRAPPERS
#define STATS_ENABLE_EIGEN_WRAPPERS
#include "stats.hpp"
#include "Eigen/Dense"
#include <iostream>

int main()
{
    Eigen::Vector<double, 9> probabilities;
    probabilities << 0.24, 0.74, 0.12, 0.56, 0.81, 0.92, 0.65, 0.57, 0.35;
    for (int32_t dof = 1; dof < 10; dof++)
        auto estimated = stats::qchisq(probabilities, static_cast<double>(dof));
    return 0;
}

Only in MSVC, template type deduction fails. Following error occurs.

error C2672: 'stats::qchisq': no matching overloaded function found
3>D:\var\stats\include\stats\stats_incl\quant\qchisq.hpp(156,1): message : could be 'Eigen::Matrix<rT,iTr,iTc,0|_Rows==1&&_Cols!=1?Eigen::RowMajor:_Cols==1&&_Rows!=1?Eigen::ColMajor:Eigen::ColMajor,_Rows,_Cols> stats::qchisq(const Eigen::Matrix<eT,iTr,iTc,0|_Rows==1&&_Cols!=1?Eigen::RowMajor:_Cols==1&&_Rows!=1?Eigen::ColMajor:Eigen::ColMajor,_Rows,_Cols> &,const T1)'
message : 'Eigen::Matrix<rT,iTr,iTc,0|_Rows==1&&_Cols!=1?Eigen::RowMajor:_Cols==1&&_Rows!=1?Eigen::ColMajor:Eigen::ColMajor,_Rows,_Cols> stats::qchisq(const Eigen::Matrix<eT,iTr,iTc,0|_Rows==1&&_Cols!=1?Eigen::RowMajor:_Cols==1&&_Rows!=1?Eigen::ColMajor:Eigen::ColMajor,_Rows,_Cols> &,const T1)': could not deduce template argument for 'const T1' from 'double'
3>D:\var\stats\include\stats\stats_incl\quant\qchisq.hpp(70,1): message : or       'std::vector<rT,std::allocator<Char>> stats::qchisq(const std::vector<_Ty,std::allocator<_Ty>> &,const T1)'
message : 'std::vector<rT,std::allocator<Char>> stats::qchisq(const std::vector<_Ty,std::allocator<_Ty>> &,const T1)': could not deduce template argument for 'const std::vector<_Ty,std::allocator<_Ty>> &' from 'Eigen::Matrix<double,9,1,0,9,1>'
3>D:\var\stats\include\stats\stats_incl\quant\qchisq.hpp(46,1): message : or       'std::conditional<std::is_integral<std::common_type<_Ty1,_Ty2>::type>::value,double,std::common_type<_Ty1,_Ty2>::type>::type stats::qchisq(const T1,const T2) noexcept'
message : Failed to specialize function template 'std::conditional<std::is_integral<std::common_type<_Ty1,_Ty2>::type>::value,double,std::common_type<_Ty1,_Ty2>::type>::type stats::qchisq(const T1,const T2) noexcept'
message : With the following template arguments:
message : 'T1=Eigen::Matrix<double,9,1,0,9,1>'
message : 'T2=double'

Only it makes above error on this environment(but when using std::vector wrapper, it is fine).

We've been using the chi-squared distribution function, and it seems to give us wrong answers. We've compared with the libgsl implementation, and it seems that for many cases, pchisq is just wrong.

For example, test_statistic = 10605, extra_parameters = 9

double p_value1 = 1.0 - gsl_cdf_chisq_P (test_statistic, (double)extra_parameters);
double p_value2 = gsl_cdf_chisq_Q (test_statistic, (double)extra_parameters);
double p_value = 1.0 - stats::pchisq(test_statistic, extra_parameters, false);

p_value is 0, but p_value1 = 1, which is the correct answer. Any ideas what's going on? Are we using it wrong?

Edit: It seems like values lower than 0.05 are often wrong.

For some inputs (apparently close to 1), qinvgamma is returning nan.
Example:

stats::qinvgamma(0.999551553841898, 2.0, 1.0)

Probably the result of overflow?