From 7b7fb09433651d82056519257a857723b9409ecb Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Fri, 29 Dec 2023 15:01:35 +0100 Subject: [PATCH] Fix wrong calculation of median and q95 in gmtbinstats (#8243) See forum post for background. In what appears to be a copy/paste error, the calculation taking an average (media if n is even) and quantile (weighted average if quantile is not exactly on a single value) then the index of the second array value was junk. Fixing this and of course q95 is always > median; well >= in the case we only have a sample of size 1. --- src/gmtbinstats.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/gmtbinstats.c b/src/gmtbinstats.c index 1323781528f..a182b93542c 100644 --- a/src/gmtbinstats.c +++ b/src/gmtbinstats.c @@ -746,7 +746,7 @@ EXTERN_MSC int GMT_gmtbinstats (void *V_API, int mode, void *args) { gmt_M_grd_loop (GMT, Grid, rowu, colu, ij) { if ((k = n_in_circle[kk])) { gmt_sort_array (GMT, zp[kk], k, GMT_FLOAT); - Grid->data[ij] = (k%2) ? zp[kk][k/2] : 0.5 * (zp[kk][(k-1)/2] + zp[kk][rowu/2]); + Grid->data[ij] = (k%2) ? zp[kk][k/2] : 0.5 * (zp[kk][(k-1)/2] + zp[kk][k/2]); gmt_M_free (GMT, zp[kk]); n++; } @@ -787,6 +787,7 @@ EXTERN_MSC int GMT_gmtbinstats (void *V_API, int mode, void *args) { case GMTBINSTATS_QUANT: /* Compute plain quantile */ gmt_M_grd_loop (GMT, Grid, rowu, colu, ij) { if ((k = n_in_circle[kk])) { + gmt_sort_array (GMT, zp[kk], k, GMT_FLOAT); Grid->data[ij] = (gmt_grdfloat) gmt_quantile_f (GMT, zp[kk], Ctrl->C.quant, k); gmt_M_free (GMT, zp[kk]); n++; @@ -882,7 +883,7 @@ EXTERN_MSC int GMT_gmtbinstats (void *V_API, int mode, void *args) { gmt_M_grd_loop (GMT, Grid, rowu, colu, ij) { if ((k = n_in_circle[kk])) { gmt_sort_array (GMT, zp[kk], k, GMT_FLOAT); - median = (k%2) ? zp[kk][k/2] : 0.5 * (zp[kk][(k-1)/2] + zp[kk][rowu/2]); + median = (k%2) ? zp[kk][k/2] : 0.5 * (zp[kk][(k-1)/2] + zp[kk][k/2]); gmt_getmad_f (GMT, zp[kk], k, median, &mad); Grid->data[ij] = mad; gmt_M_free (GMT, zp[kk]);