/usr/lib/python2.7/dist-packages/astroML/lumfunc.py is in python-astroml 0.3-7.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 | import numpy as np
def _sorted_interpolate(x, y, x_eval):
"""utility function for binned_Cminus"""
# note that x should be sorted
N = len(x)
ind = x.searchsorted(x_eval)
ind[ind == N] = N - 1
y_eval = np.zeros(x_eval.shape)
# find perfect matches
match = (x[ind] == x_eval) | (x_eval > x[-1]) | (x_eval < x[0])
y_eval[match] = y[ind[match]]
ind = ind[~match]
# take care of extrapolation
ind[ind == 0] = 1
x_lo = x[ind - 1]
x_up = x[ind]
y_lo = y[ind - 1]
y_up = y[ind]
# take care of places where x_lo = x_up
y_eval[~match] = (y_lo + (x_eval[~match] - x_lo)
* (y_up - y_lo) / (x_up - x_lo))
return y_eval
def Cminus(x, y, xmax, ymax):
"""Lynden-Bell's C-minus method
Parameters
----------
x : array_like
array of x values
y : array_like
array of y values
xmax : array_like
array of maximum x values for each y value
ymax : array_like
array of maximum y values for each x value
Returns
-------
Nx, Ny, cuml_x, cuml_y: ndarrays
Nx and cuml_x are in the order of the sorted x array
Ny and cuml_y are in the order of the sorted y array
"""
# make copies of input
x, y, xmax, ymax = map(np.array, (x, y, xmax, ymax))
Nall = len(x)
cuml_x = np.zeros(x.shape)
cuml_y = np.zeros(y.shape)
Nx = np.zeros(x.shape)
Ny = np.zeros(y.shape)
# first the y direction.
i_sort = np.argsort(y)
x = x[i_sort]
y = y[i_sort]
xmax = xmax[i_sort]
ymax = ymax[i_sort]
for j in range(1, Nall):
Ny[j] = np.sum(x[:j] < xmax[j])
Ny[0] = np.inf
cuml_y = np.cumprod(1. + 1. / Ny)
Ny[0] = 0
# renormalize
cuml_y *= Nall / cuml_y[-1]
#now the x direction
i_sort = np.argsort(x)
x = x[i_sort]
y = y[i_sort]
xmax = xmax[i_sort]
ymax = ymax[i_sort]
for i in range(1, Nall):
Nx[i] = np.sum(y[:i] < ymax[i])
Nx[0] = np.inf
cuml_x = np.cumprod(1. + 1. / Nx)
Nx[0] = 0
# renormalize
cuml_x *= Nall / cuml_x[-1]
return Nx, Ny, cuml_x, cuml_y
def binned_Cminus(x, y, xmax, ymax, xbins, ybins, normalize=False):
"""Compute the binned distributions using the Cminus method
Parameters
----------
x : array_like
array of x values
y : array_like
array of y values
xmax : array_like
array of maximum x values for each y value
ymax : array_like
array of maximum y values for each x value
xbins : array_like
array of bin edges for the x function: size=Nbins_x + 1
ybins : array_like
array of bin edges for the y function: size=Nbins_y + 1
normalize : boolean
if true, then returned distributions are normalized. Default
is False.
Returns
-------
dist_x, dist_y : ndarrays
distributions of size Nbins_x and Nbins_y
"""
Nx, Ny, cuml_x, cuml_y = Cminus(x, y, xmax, ymax)
# simple linear interpolation using a binary search
# interpolate the cumulative distributions
x_sort = np.sort(x)
y_sort = np.sort(y)
Ix_edges = _sorted_interpolate(x_sort, cuml_x, xbins)
Iy_edges = _sorted_interpolate(y_sort, cuml_y, ybins)
if xbins[0] < x_sort[0]:
Ix_edges[0] = cuml_x[0]
if xbins[-1] > x_sort[-1]:
Ix_edges[-1] = cuml_x[-1]
if ybins[0] < y_sort[0]:
Iy_edges[0] = cuml_y[0]
if ybins[-1] > y_sort[-1]:
Iy_edges[-1] = cuml_y[-1]
x_dist = np.diff(Ix_edges) / np.diff(xbins)
y_dist = np.diff(Iy_edges) / np.diff(ybins)
if normalize:
x_dist /= len(x)
y_dist /= len(y)
return x_dist, y_dist
def bootstrap_Cminus(x, y, xmax, ymax, xbins, ybins,
Nbootstraps=10, normalize=False):
"""
Compute the binned distributions using the Cminus method, with
bootstrapped estimates of the errors
Parameters
----------
x : array_like
array of x values
y : array_like
array of y values
xmax : array_like
array of maximum x values for each y value
ymax : array_like
array of maximum y values for each x value
xbins : array_like
array of bin edges for the x function: size=Nbins_x + 1
ybins : array_like
array of bin edges for the y function: size=Nbins_y + 1
Nbootstraps : int
number of bootstrap resamplings to perform
normalize : boolean
if true, then returned distributions are normalized. Default
is False.
Returns
-------
dist_x, err_x, dist_y, err_y : ndarrays
distributions of size Nbins_x and Nbins_y
"""
x, y, xmax, ymax = map(np.asarray, (x, y, xmax, ymax))
x_dist = np.zeros((Nbootstraps, len(xbins) - 1))
y_dist = np.zeros((Nbootstraps, len(ybins) - 1))
for i in range(Nbootstraps):
ind = np.random.randint(0, len(x), len(x))
x_dist[i], y_dist[i] = binned_Cminus(x[ind], y[ind],
xmax[ind], ymax[ind],
xbins, ybins,
normalize=normalize)
return (x_dist.mean(0), x_dist.std(0, ddof=1),
y_dist.mean(0), y_dist.std(0, ddof=1))
|