# [SciPy-User] scipy.stats.nanmedian

josef.pktd@gmai... josef.pktd@gmai...
Fri Jan 22 15:08:23 CST 2010

```On Fri, Jan 22, 2010 at 12:03 PM,  <josef.pktd@gmail.com> wrote:
> On Fri, Jan 22, 2010 at 11:52 AM, Keith Goodman <kwgoodman@gmail.com> wrote:
>> On Fri, Jan 22, 2010 at 8:46 AM,  <josef.pktd@gmail.com> wrote:
>>> On Fri, Jan 22, 2010 at 11:09 AM, Keith Goodman <kwgoodman@gmail.com> wrote:
>>>> On Thu, Jan 21, 2010 at 8:18 PM,  <josef.pktd@gmail.com> wrote:
>>>>> On Thu, Jan 21, 2010 at 10:01 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
>>>>>> On Thu, Jan 21, 2010 at 6:41 PM, Pierre GM <pgmdevlist@gmail.com> wrote:
>>>>>>> On Jan 21, 2010, at 9:28 PM, Keith Goodman wrote:
>>>>>>>> That's the only was I was able to figure out how to pull 1.0 out of
>>>>>>>> np.array(1.0). Is there a better way?
>>>>>>>
>>>>>>>
>>>>>>> .item()
>>>>>>
>>>>>> Thanks. item() looks better than tolist().
>>>>>>
>>>>>> I simplified the function:
>>>>>>
>>>>>> def nanmedian(x, axis=0):
>>>>>>    x, axis = _chk_asarray(x,axis)
>>>>>>    if x.ndim == 0:
>>>>>>        return float(x.item())
>>>>>>    x = x.copy()
>>>>>>    x = np.apply_along_axis(_nanmedian,axis,x)
>>>>>>    if x.ndim == 0:
>>>>>>        x = float(x.item())
>>>>>>    return x
>>>>>>
>>>>>> and opened a ticket:
>>>>>>
>>>>>> http://projects.scipy.org/scipy/ticket/1098
>>>>>
>>>>>
>>>>> How about getting rid of apply_along_axis?    see attachment
>>>>>
>>>>> I don't know whether or how much faster it is, but there is a ticket
>>>>> that the current version is slow.
>>>>> No hidden bug or corner case guarantee yet.
>>>>
>>>> It is faster. But here is one case it does not handle:
>>>>
>>>>>> nanmedian([1, 2])
>>>>   array([ 1.5])
>>>>>> np.median([1, 2])
>>>>   1.5
>>>>
>>>> I'm sure it could be fixed. But having to fix it, and the fact that it
>>>> is a larger change, decreases the likelihood that it will make it into
>>>> the next version of scipy. One option is to make the small bug fix I
>>>> suggested (ticket #1098) and add the corresponding unit tests. Then we
>>>> can take our time to design a better version of nanmedian.
>>>
>>> I didn't see the difference to np.median for this case, I think I was
>>> taking the shape answer from the other thread on the return of splines
>>> and interpolation.
>>>
>>> If I change the last 3 lines to
>>>    if nanmed.size == 1:
>>>       return nanmed.item()
>>>    return nanmed
>>>
>>> then I get agreement with numpy for the following test cases
>>>
>>> print nanmedian(1), np.median(1)
>>> print nanmedian(np.array(1)), np.median(1)
>>> print nanmedian(np.array([1])), np.median(np.array([1]))
>>> print nanmedian(np.array([[1]])), np.median(np.array([[1]]))
>>> print nanmedian(np.array([1,2])), np.median(np.array([1,2]))
>>> print nanmedian(np.array([[1,2]])), np.median(np.array([[1,2]]),axis=0)
>>> print nanmedian([1]), np.median([1])
>>> print nanmedian([[1]]), np.median([[1]])
>>> print nanmedian([1,2]), np.median([1,2])
>>> print nanmedian([[1,2]]), np.median([[1,2]],axis=0)
>>> print nanmedian([1j,2]), np.median([1j,2])
>>>
>>> Am I still missing any cases?
>>>
>>> The vectorized version should be faster for this case
>>> http://projects.scipy.org/scipy/ticket/740
>>> but maybe not for long and narrow arrays.
>>
>> Here is an odd one:
>>
>>>> nanmedian(True)
>>   1.0
>>>> nanmedian([True])
>>   0.5  # <--- strange
>>
>>>> np.median(True)
>>   1.0
>>>> np.median([True])
>>   1.0
>
> definitely weird
>
>>>> (np.array(True)+np.array(True))/2.
> 0.5
>>>> np.array([True, True]).sum()
> 2
>>>> np.array([True, True]).mean()
> 1.0
>
> I assumed mean (is used by np.ma.median) is the same as adding and dividing by 2
>
> Josef
>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>

It got a bit ugly, too many shapes for a single number, and tricky
axis handling.
Your idea of making just the smallish changes looks more attractive now.

and almost correct

>>> np.median(np.array([[[1]]]),axis=0).shape
(1, 1)
>>> nanmedian(np.array([[[1]]]),axis=0).shape
(1, 1)
>>> np.median(np.array([[[1]]]),axis=-1).shape
(1, 1)
>>> nanmedian(np.array([[[1]]]),axis=-1).shape
(1, 1)

but there slipped a python in

>>> nanmedian(np.array([[[1]]]),axis=None).__class__
<type 'float'>
>>> np.median(np.array([[[1]]]),axis=None).__class__
<type 'numpy.float64'>

.item() returns a python number not a numpy number

>>> np.array([[[1]]]).item().__class__
<type 'int'>
>>> np.array([[[1]]]).flat[0].__class__
<type 'numpy.int32'>

I also didn't know this:

>>> None < 0
True

Josef
-------------- next part --------------
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 20 10:18:32 2010
Author: josef-pktd
"""

import numpy as np
from scipy import stats

def nanmedian(x, axis = 0):
keepshape = list(np.shape(x))
x, axis2 = stats.stats._chk_asarray(x, axis)
if (not axis is None) and axis<0 : # and x.ndim>2:
axis = x.ndim + axis
#print 'axis', axis

#print x, keepshape
if x.ndim == 0 or (x.size==1 and axis is None):
return 1.0*x.item()
if keepshape and not axis is None :
keepshape.pop(axis)
if x.size == 1: #x.ndim == 0:
return 1.0*x.reshape(keepshape)
if x.dtype == np.bool:
x = x.astype(int)
if axis is None:
axis = 0

x = np.sort(x, axis=axis)
nall = x.shape[axis]
notnancount = nall - np.isnan(x).sum(axis=axis)
#print 'notnancount', notnancount

(idx, rmd) = divmod(notnancount, 2)
#idx = np.atleast_1d(idx)
#print 'idx', idx
#print idx.shape
indx = np.ogrid[[slice(i)for i in x.shape]]
indxlo = indx[:]
idxslice = map(slice, idx.shape)
idxslice.insert(axis, None)
#print idxslice
idx = np.atleast_1d(idx)
indxlo[axis] = idx[idxslice]
#print indxlo
indxhi = indx[:]
indxhi[axis] = (idx - (1-rmd))[idxslice]#[idx.shape[:axis]+(None,)+idx.shape[axis:]]
#print map(np.shape, indxhi)
#print indxhi
nanmed = (x[indxlo] + x[indxhi])/2.

#print 'keepshape, nanmed.shape',keepshape, nanmed.shape
if np.ndim == 0:
return nanmed.item()
#return nanmed.reshape(keepshape)
return np.squeeze(nanmed) #.reshape(keepshape)
if nanmed.size == 1:
return nanmed.reshape(keepshape)
return nanmed.item()
return nanmed

from numpy.testing import assert_equal, assert_almost_equal

for axis in [0,1, None, -1]:
for i in range(5):
# for complex
#x = 1j+np.arange(20).reshape(4,5)
x = np.arange(20).reshape(4,5).astype(float)
x[zip(np.random.randint(4, size=(2,5)))] = np.nan
assert_equal(nanmedian(x, axis=0), stats.nanmedian(x, axis=0))

for axis in [0,1,2, None, -1]:
for i in range(5):
x = np.arange(3*4*5).reshape(3,4,5).astype(float)
x[np.random.randint(3, size=(3,4,5))] = np.nan
assert_equal(nanmedian(x, axis=0), stats.nanmedian(x, axis=0))

xli = [ [1], [[1]], [1,2], [1j], [1j,2], [True], [False],
[True,False], [True, False, True], np.round(np.random.randn(2,4,5),4)]
xxli = xli + map(np.array,xli)

for axis in [0, -1, None]:
print '\n next case',axis
for x in xxli:
try:
assert_equal(nanmedian(x, axis=axis), np.median(x, axis=axis))
assert_equal(np.shape(nanmedian(x, axis=axis)), np.shape(np.median(x, axis=axis)))
except:
print 'failure with', x
print nanmedian(x, axis=axis),  np.median(x, axis=axis)
raise

y=np.round(np.random.randn(2,3,5),4)
axis = -1 #None
print np.median(y,axis=axis)
nm = nanmedian(y,axis=axis)
print nm
print np.median(y,axis=axis) - nm

#np.broadcast(array([[[0]],[[1]]]), array([[[1, 1, 1, 1, 1]], [[1, 1, 1, 1, 1]]]), array([[[0, 1, 2, 3, 4]]]))

```