From 20dd57064a70b3702b04043498c9b583ee130b43 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Fri, 29 Dec 2023 14:39:14 +0100 Subject: [PATCH] Fix wrong calculation of median and q95 in gmtbinstats 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]);