1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
import numpy as N
import numpy.linalg as L
import scipy.interpolate
import scipy.linalg
__docformat__ = 'restructuredtext'
def recipr(X):
"""
Return the reciprocal of an array, setting all entries less than or
equal to 0 to 0. Therefore, it presumes that X should be positive in
general.
"""
x = N.maximum(N.asarray(X).astype(N.float64), 0)
return N.greater(x, 0.) / (x + N.less_equal(x, 0.))
def mad(a, c=0.6745, axis=0):
"""
Median Absolute Deviation:
median(abs(a - median(a))) / c
"""
_shape = a.shape
a.shape = N.product(a.shape,axis=0)
m = N.median(N.fabs(a - N.median(a))) / c
a.shape = _shape
return m
def recipr0(X):
"""
Return the reciprocal of an array, setting all entries equal to 0
as 0. It does not assume that X should be positive in
general.
"""
test = N.equal(N.asarray(X), 0)
return N.where(test, 0, 1. / X)
def clean0(matrix):
"""
Erase columns of zeros: can save some time in pseudoinverse.
"""
colsum = N.add.reduce(matrix**2, 0)
val = [matrix[:,i] for i in N.flatnonzero(colsum)]
return N.array(N.transpose(val))
def rank(X, cond=1.0e-12):
"""
Return the rank of a matrix X based on its generalized inverse,
not the SVD.
"""
D = scipy.linalg.svdvals(X)
return int(N.add.reduce(N.greater(D / D.max(), cond).astype(N.int32)))
def fullrank(X, r=None):
"""
Return a matrix whose column span is the same as X.
If the rank of X is known it can be specified as r -- no check
is made to ensure that this really is the rank of X.
"""
if r is None:
r = rank(X)
V, D, U = L.svd(X, full_matrices=0)
order = N.argsort(D)
order = order[::-1]
value = []
for i in range(r):
value.append(V[:,order[i]])
return N.asarray(N.transpose(value)).astype(N.float64)
class StepFunction:
"""
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
--------
>>> from numpy import arange
>>> from scipy.sandbox.models.utils import StepFunction
>>>
>>> x = arange(20)
>>> y = arange(20)
>>> f = StepFunction(x, y)
>>>
>>> print f(3.2)
3.0
>>> print f([[3.2,4.5],[24,-3.1]])
[[ 3. 4.]
[ 19. 0.]]
"""
def __init__(self, x, y, ival=0., sorted=False):
_x = N.asarray(x)
_y = N.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 = N.hstack([[-N.inf], _x])
self.y = N.hstack([[ival], _y])
if not sorted:
asort = N.argsort(self.x)
self.x = N.take(self.x, asort, 0)
self.y = N.take(self.y, asort, 0)
self.n = self.x.shape[0]
def __call__(self, time):
tind = N.searchsorted(self.x, time) - 1
_shape = tind.shape
return self.y[tind]
def ECDF(values):
"""
Return the ECDF of an array as a step function.
"""
x = N.array(values, copy=True)
x.sort()
x.shape = N.product(x.shape,axis=0)
n = x.shape[0]
y = (N.arange(n) + 1.) / n
return StepFunction(x, y)
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 = N.array(y)
a = N.argsort(y)
return scipy.interpolate.interp1d(y[a], x[a])
|
