Skip to content

wrong result with xt::quantile #2651

Description

@ThibHlln

Hi,

I seem to be getting the wrong result with xt::quantile in some occasions (i.e. with some specific data like the one below, and only on linux with GCC so far).

For example:

I need to pre-compute my quantiles based on some confidence levels:

// pre-compute quantiles based on confidence levels
xt::xtensor<double, 1> lvl = {30., 80.};
xt::xtensor<double, 2> qth = xt::zeros<double>({n_itv, std::size_t {2}});
xt::col(qth, 0) = 0.5 - lvl / 200.; 
xt::col(qth, 1) = 0.5 + lvl / 200.; 
std::cout << qth << std::endl;
// {{0.35, 0.1 },
//  {0.65, 0.9 }};

// "hand-written" quantiles
xt::xtensor<double, 2> qth_ = {{0.35, 0.10},
                               {0.65, 0.90}};

Then, when I use these pre-computed quantiles with xt::quantile I don't get the right result:

// on 2D tensors
xt::xtensor<double, 2> prd = {{927., 932., 935., 943., 944., 944., 945., 948., 952.,
                               962., 968., 968., 969., 969., 974., 981., 993., 994.,
                               994., 1003., 1007., 1008., 1008., 1012., 1013., 1014., 1014.,
                               1017., 1019., 1020., 1023., 1026., 1026., 1028., 1030., 1032.,
                               1039., 1047., 1071., 1080., 1085., 1088., 1111., 1112., 1117.,
                               1119., 1128., 1130., 1209., 1309., 1426.}};

// wrong result if I pass a view of pre-computed quantiles
std::cout << xt::quantile(prd, xt::view(qth, 0), 1) << std::endl;
// {{ 1000.5},
//  { 1027. }}

// right result if I pass quantiles directly in function call
std::cout << xt::quantile(prd, {0.35, 0.1}, 1) << std::endl;
// {{ 994.},
//  { 944.}}

// right result if I pass a view of "hand-written" quantiles
std::cout << xt::quantile(prd, xt::view(qth_, 0), 1) << std::endl;
// {{ 994.},
//  { 944.}}

And if I work on 1D tensors instead, I get a result half wrong:

// on 1D tensors
xt::xtensor<double, 1> prd = {927., 932., 935., 943., 944., 944., 945., 948., 952.,
                              962., 968., 968., 969., 969., 974., 981., 993., 994.,
                              994., 1003., 1007., 1008., 1008., 1012., 1013., 1014., 1014.,
                              1017., 1019., 1020., 1023., 1026., 1026., 1028., 1030., 1032.,
                              1039., 1047., 1071., 1080., 1085., 1088., 1111., 1112., 1117.,
                              1119., 1128., 1130., 1209., 1309., 1426.};

// wrong result for 2nd quantile if I pass a view of pre-computed quantiles
std::cout << xt::quantile(prd, xt::view(qth, 0))  << std::endl;
// {  994.,  1027.}

// right result if I pass quantiles directly in function call
std::cout << xt::quantile(prd {0.35, 0.1}) << std::endl;
// { 994.,  944.}

// right result if I pass a view of "hand-written" quantiles
std::cout << xt::quantile(prd, xt::view(qth_, 0)) << std::endl;
// { 994.,  944.}

I am unsure whether I am doing something wrong or if this is a bug?

I first thought it was because I was passing views for the probas parameter of xt::quantile, but as shown above this works with the "hand-written" tensor.

I also tried adding xt::eval when pre-computing my quantiles, but it does not seem to solve the problem.

I would appreciate any help/ideas on this problem.

Thanks.


For reference, with numpy:

prd = numpy.array([[927.,  932.,  935.,  943.,  944.,  944.,  945.,  948.,  952.,
                    962.,  968.,  968.,  969.,  969.,  974.,  981.,  993.,  994.,
                    994., 1003., 1007., 1008., 1008., 1012., 1013., 1014., 1014.,
                    1017., 1019., 1020., 1023., 1026., 1026., 1028., 1030., 1032.,
                    1039., 1047., 1071., 1080., 1085., 1088., 1111., 1112., 1117.,
                    1119., 1128., 1130., 1209., 1309., 1426.]])

print(numpy.quantile(prd, [0.35, 0.1], 1))
# [[994.]
#  [944.]]

Activity

  1. AntoinePrv commented on Feb 8, 2023

    @AntoinePrv
    Contributor

    Hello @ThibHlln , thanks for reporting!

    Looking into your 1st example, if you increase the cout precision before your example with std::cout.precision(17);, I notice that:

    std::cout << qth << std::endl;
    // {0.34999999999999998, 0.65000000000000002},
    // {0.09999999999999998, 0.90000000000000002}}

    While qth_ still prints the same

    std::cout << qth_ << std::endl;
    // {{ 0.35,  0.1 },
    // { 0.65,  0.9 }}

    If I use qth_ with these values

    // "hand-written" quantiles
    xt::xtensor<double, 2> qth_ = {{0.34999999999999998, 0.65000000000000002},,
                                   {0.09999999999999998, 0.90000000000000002}};

    I get the same result {{ 1000.5}, { 1027. }} in all cases.

    Now the weird part, is at the double precision, 0.34999999999999998 and 0.35 are the same number! Same goes for the others.
    The first question for me is why qth_ prints differently depending on whether I pass it 0.34999999999999998 or 0.35?

  2. AntoinePrv commented on Feb 8, 2023

    @AntoinePrv
    Contributor

    @ThibHlln forget everything I said.

    qth is not what you say it is, it is the transposed of qth_.

    std::cout << qth << std::endl;
    // {{0.35, 0.65 },
    //  {0.1, 0.9 }};
    

    The "weird" part I was seeing was an artifact of copy-pasting without noticing one was the transposed of the other.

  3. ThibHlln commented on Feb 9, 2023

    @ThibHlln
    ContributorAuthor

    Hi @AntoinePrv ,

    Thanks a lot for looking into this. Apologies for the mistake in my MWE.
    Nevertheless, at location [0, 0] we still have 0.35, right? So the answer should not be 1000.5 but 994?

  4. AntoinePrv commented on Feb 9, 2023

    @AntoinePrv
    Contributor

    If I'm not wrong, we do have the consistent results regardless of whether we use a view, or an hand-written array.

    Nevertheless, at location [0, 0] we still have 0.35, right? So the answer should not be 1000.5 but 994?

    Indeed. So it looks like like that the function gives a different output depending on whether the proba afterwards is larger or smaller, which is definitely a bug! I'll investigate.

  5. ThibHlln commented on Feb 9, 2023

    @ThibHlln
    ContributorAuthor

    If I'm not wrong, we do have the consistent results regardless of whether we use a view, or an hand-written array.

    Yes, absolutely. My diagnosis was wrong.

    Indeed. So it looks like like that the function gives a different output depending on whether the proba afterwards is larger or smaller, which is definitely a bug! I'll investigate.

    Thank you very much!

  6. AntoinePrv commented on Feb 9, 2023

    @AntoinePrv
    Contributor

    @ThibHlln I think I am getting starting to get confused with the example... Right now I don't get any wrong output.
    Can you do a fresh pass on it and confirm what answer is bug you are getting.

  7. ThibHlln commented on Feb 9, 2023

    @ThibHlln
    ContributorAuthor

    Sure, here is an updated MWE:

    xt::xtensor<double, 2> prd = {{927., 932., 935., 943., 944., 944., 945., 948., 952.,
                                   962., 968., 968., 969., 969., 974., 981., 993., 994.,
                                   994., 1003., 1007., 1008., 1008., 1012., 1013., 1014., 1014.,
                                   1017., 1019., 1020., 1023., 1026., 1026., 1028., 1030., 1032.,
                                   1039., 1047., 1071., 1080., 1085., 1088., 1111., 1112., 1117.,
                                   1119., 1128., 1130., 1209., 1309., 1426.}};
    
    std::cout << xt::quantile(prd, {0.10}, 1) << std::endl;
    // {{ 944.}}
    
    std::cout << xt::quantile(prd, {0.35}, 1) << std::endl;
    // {{ 994.}}
    
    std::cout << xt::quantile(prd, {0.65}, 1) << std::endl;
    // {{ 1027.}}
    
    std::cout << xt::quantile(prd, {0.35, 0.10}, 1) << std::endl;
    // {{ 994.},
    //  { 944.}}
    
    // WRONG RESULT FOR PROBA=0.35
    std::cout << xt::quantile(prd, {0.35, 0.65}, 1) << std::endl;
    // {{ 1000.5},  <<<< ?!?
    //  { 1027. }}
    
    
    // WRONG RESULT FOR PROBA=0.35
    std::cout << xt::quantile(prd, {0.65, 0.35}, 1) << std::endl;
    // {{ 1027. }, 
    //  { 1000.5}} <<<< ?!?
    
  8. AntoinePrv commented on Feb 10, 2023

    @AntoinePrv
    Contributor

    Thanks @ThibHlln. The issue seems to come from partition, for which I opened a PR.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions