Describe the issue:

nanmedian will sometimes return inf instead of a proper value for numbers that are close the maximum value of floats data type.
This behaviour seems to depend on the usage of the "axis" keyword parameter.

For example np.nanmedian(np.array([[1,4e4]],dtype="float16"),axis=0)
returns [ 1., inf] instead of [1.,4e4]

This is true at least for float16, float32 and float64 with values that are close to the respective max of these data types.

np.median doesn't seem have any issues with these values.

Reproduce the code example:

import numpy as np

# I'll start with a simple array
# max value for float16 should be 6.55e4
a = np.array([1,1e4,4e4],dtype="float16")
print(a)

# then try to compute the median along a given axis
# this is the line that produces faulty output
b = np.nanmedian([a,a],axis=0)
print(b)
# we get [1.e+00 1.e+04    inf] instead of [1.e+00 1.e+04 4.e+04]


# the regular median shows the expected behaviour
c = np.median([a,a],axis=0)
print(c)


# the strange behaviour is not always present
# for example nanmedian works ok in this situation
d = np.nanmedian(np.array([[4e4,4e4]],dtype="float16"))
print(d)

Error message:

No error message but a runtime warning :
C:\Users\***\Anaconda3\lib\site-packages\numpy\core\_methods.py:48: RuntimeWarning: overflow encountered in reduce
  return umr_sum(a, axis, dtype, out, keepdims, initial, where)

NumPy/Python version information:

1.21.5 3.9.12 (main, Apr 4 2022, 05:22:27) [MSC v.1916 64 bit (AMD64)]

Context for the issue:

This issue causes valid values that are not inf to be treated as inf (and eventually as NaNs after further operations) unless a larger data type is used which isn't always feasible when treating very large data sets.