wrong result with xt::quantile

[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.]]

[AntoinePrv]

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?

[AntoinePrv]

@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.

[ThibHlln]

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?

[AntoinePrv]

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.

[ThibHlln]

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!

[AntoinePrv]

@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.

[ThibHlln]

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}} <<<< ?!?

[AntoinePrv]

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