Skip to content

Commit

Permalink
Fix wrong calculation of median and q95 in gmtbinstats (#8243)
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
PaulWessel authored Dec 29, 2023
1 parent 3199403 commit 7b7fb09
Showing 1 changed file with 3 additions and 2 deletions.
5 changes: 3 additions & 2 deletions src/gmtbinstats.c
Original file line number Diff line number Diff line change
Expand Up @@ -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++;
}
Expand Down Expand Up @@ -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++;
Expand Down Expand Up @@ -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]);
Expand Down

0 comments on commit 7b7fb09

Please sign in to comment.