search

pacificclimate / climdex.pcic

Routines to compute ETCCDI Climdex indices
GNU General Public License v3.0
23 stars 13 forks source link

c_quantile may use the wrong formula for nppm #24

Closed bzah closed 3 years ago

bzah commented 3 years ago

I'm working on a similar libraries in python (Icclim/Xclim) and I'm currently on the computation of percentiles. I noticed you rewrote the 8th method of Hyndman&Fan from R but I found an unusual -1 in your formula so I wanted to notice you. In _zhang_running_quantile.cc::cquantile method. you compute the virtual index of the percentile with const double nppm = a + quantile * (n + 1 - a - b) - 1; Because a and b are always 1/3 in this case this can be simplified:

nppm = 1/3 + quantile * (n + 1 - 1/3 - 1/3) - 1;
nppm = quantile * (n + 1/3) - 2/3;

This is quite different from the other implementation I came across (Scipy, R and Icclim in it's current state). When developed, they all came down to:

nppm = quantile * (n + 1/3) + 1/3

If this behavior is not intended, I would suggest to simply remove the -1 in the formula.

For reference you can find theses implementations there: Scipy R Icclim

bzah commented 3 years ago

Well R arrays start at 1 but in most other languages it start at 0. -1 is here exactly for that. So it seems Climdex is right and Scipy and Icclim are not.