quantile_curve(curves, quantile, weights) is broken when weights is None about oq-hazardlib HOT 9 CLOSED

I agree we must be consistent in the way we calculate quantiles hence I would eliminate the 'customised' implementation

from oq-hazardlib.

When the individual curves are obtained by sampling the logic-tree and no weights are assigned, we do not know the empirical distribution of the sampled values. To estimate statistics from the sampled values, we would need to assign a plotting position based on assumptions regarding the underlying distribution and depending upon the statistic we are trying to estimate. The plotting position recommended by Cunnae (1978) for unbiased estimation of quantiles when the underlying distribution type is unknown seems to be the one implemented in the snippet above. The equation recommended in that paper by Cunnae is listed here.

We don't have to go through the same process to estimate the mean because the sample mean is already an unbiased estimator.

from oq-hazardlib.

Here is an example that uses the algorithm currently implemented for the full enumeration case:

  quantile = 0.75
  poes = numpy.array([.99, .95, .93, .87, .60])
  weights = numpy.array([.4, .1, .3, .1, .1])
  sorted_poe_idxs = numpy.argsort(poes)  # array([4, 3, 2, 1, 0])
  sorted_poes = poes[sorted_poe_idxs]  # array([ 0.6 ,  0.87,  0.93,  0.95,  0.99])
  sorted_weights = weights[sorted_poe_idxs]  # array([ 0.1,  0.1,  0.3,  0.1,  0.4])
  cum_weights = numpy.cumsum(sorted_weights)  # array([ 0.1,  0.2,  0.5,  0.6,  1. ])
  numpy.interp(quantile, cum_weights, sorted_poes)  # 0.965

Does anybody have an idea of the origin of this algorithm? It would be nice to have a link to some standard reference.

from oq-hazardlib.

I believe that's obtaining the quantiles simply by linear interpolation on the empirical CDF of the poes. For reference, see the Type 4 quantile as computed by R.

from oq-hazardlib.

Yes, it makes sense since it does not make any assumption on the underlying distribution. I would keep this algorithm, since it is compatible with the past. BTW, even if we remove the other algorithm for weights=None there is no impact on the risk tests, since there are no risk tests for sampling.

from oq-hazardlib.

Thanks this is all very intesting. I have (at least) one question. Isn't the scipy function taking into account the fact case where data doesn't have weight associated?

from oq-hazardlib.

It is used in that case. The issue is that the scipy function gives different numbers than the interpolated CDF approach when the weights are all equal. That make the calculation of quantiles surprising and inconsistent with the mean computation.

from oq-hazardlib.

I'll prepare some additional tests on sampling. This is also interesting: https://phobson.github.io/mpl-probscale/tutorial/closer_look_at_plot_pos.html

from oq-hazardlib.

If you want to see the effect on the current tests, see gem/oq-engine#2492

from oq-hazardlib.