# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
''' General matrix and other utilities for statistics '''
from __future__ import print_function
from __future__ import absolute_import
__docformat__ = 'restructuredtext'
import numpy as np
import scipy.interpolate
[docs]def mad(a, c=0.6745, axis=0):
"""
Median Absolute Deviation:
median(abs(a - median(a))) / c
"""
_shape = a.shape
a.shape = np.product(a.shape, axis=0)
m = np.median(np.fabs(a - np.median(a))) / c
a.shape = _shape
return m
[docs]class StepFunction(object):
""" A basic step function
Values at the ends are handled in the simplest way possible: everything to
the left of ``x[0]`` is set to `ival`; everything to the right of ``x[-1]``
is set to ``y[-1]``.
Examples
--------
>>> x = np.arange(20)
>>> y = np.arange(20)
>>> f = StepFunction(x, y)
>>>
>>> f(3.2)
3.0
>>> res = f([[3.2, 4.5],[24, -3.1]])
>>> np.all(res == [[ 3, 4],
... [19, 0]])
True
"""
[docs] def __init__(self, x, y, ival=0., sorted=False):
_x = np.asarray(x)
_y = np.asarray(y)
if _x.shape != _y.shape:
raise ValueError(
'in StepFunction: x and y do not have the same shape')
if len(_x.shape) != 1:
raise ValueError('in StepFunction: x and y must be 1-dimensional')
self.x = np.hstack([[- np.inf], _x])
self.y = np.hstack([[ival], _y])
if not sorted:
asort = np.argsort(self.x)
self.x = np.take(self.x, asort, 0)
self.y = np.take(self.y, asort, 0)
self.n = self.x.shape[0]
def __call__(self, time):
tind = np.searchsorted(self.x, time) - 1
return self.y[tind]
[docs]def ECDF(values):
"""
Return the ECDF of an array as a step function.
"""
x = np.array(values, copy=True)
x.sort()
x.shape = np.product(x.shape, axis=0)
n = x.shape[0]
y = (np.arange(n) + 1.) / n
return StepFunction(x, y)
[docs]def monotone_fn_inverter(fn, x, vectorized=True, **keywords):
"""
Given a monotone function x (no checking is done to verify monotonicity)
and a set of x values, return an linearly interpolated approximation
to its inverse from its values on x.
"""
if vectorized:
y = fn(x, **keywords)
else:
y = []
for _x in x:
y.append(fn(_x, **keywords))
y = np.array(y)
a = np.argsort(y)
return scipy.interpolate.interp1d(y[a], x[a])