import numpy as np
import scipy as sp
I. Integration and interpolation
zz = np.loadtxt('wiggleZ_DR1_z.dat',dtype='float'); # Load WiggleZ redshifts
np.min(zz)
1.0000000000000001e-05
np.max(zz)
1.9903999999999999
nbins = 50;
n, bins, patches = hist(zz,nbins)
x = bins[0:nbins] + (bins[2]-bins[1])/2; # Adjust bin edges to centre
from scipy.interpolate import interp1d
zdist = interp1d(x,n, kind='cubic', bounds_error = False, fill_value=0)
z = linspace(0,2,100)
plot(z,zdist(z))
[<matplotlib.lines.Line2D at 0x1056bfb10>]
from scipy import integrate
cumz = lambda z0: integrate.quad(zdist,0,z0)[0]
total = cumz(5)
/usr/local/Cellar/python/2.7.3/lib/python2.7/site-packages/scipy/integrate/quadpack.py:288: UserWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. warnings.warn(msg)
total
3232.7938146343804
Pz = 0 * z;
for i in range(len(z)):
Pz[i] = cumz(z[i]) / total
plot(z,Pz,Pz,z)
[<matplotlib.lines.Line2D at 0x10589a810>, <matplotlib.lines.Line2D at 0x10589aa10>]
zvar = interp1d(Pz,z,kind='cubic')
from numpy import random
rand()
0.6928841274449873
output = zvar(rand(100000,1))
hist(output,100)
(array([ 53, 28, 14, 11, 15, 10, 9, 10, 13, 6, 14, 5, 194, 129, 114, 419, 955, 922, 944, 1097, 1253, 1376, 1414, 1726, 1866, 2003, 2384, 2661, 2687, 2960, 3473, 3684, 3535, 3703, 4113, 4146, 4151, 4325, 4348, 4017, 3586, 3273, 2906, 2747, 2686, 2376, 2136, 1935, 1711, 1646, 1497, 1233, 924, 786, 787, 730, 607, 452, 363, 371, 277, 268, 248, 184, 173, 145, 66, 69, 93, 89, 66, 29, 16, 25, 35, 42, 27, 28, 33, 42, 55, 27, 35, 35, 58, 33, 30, 12, 23, 21, 22, 42, 18, 20, 32, 17, 12, 11, 2, 1]), array([-0.23761697, -0.21519182, -0.19276667, -0.17034152, -0.14791637, -0.12549122, -0.10306607, -0.08064092, -0.05821577, -0.03579063, -0.01336548, 0.00905967, 0.03148482, 0.05390997, 0.07633512, 0.09876027, 0.12118542, 0.14361057, 0.16603572, 0.18846087, 0.21088601, 0.23331116, 0.25573631, 0.27816146, 0.30058661, 0.32301176, 0.34543691, 0.36786206, 0.39028721, 0.41271236, 0.43513751, 0.45756265, 0.4799878 , 0.50241295, 0.5248381 , 0.54726325, 0.5696884 , 0.59211355, 0.6145387 , 0.63696385, 0.659389 , 0.68181415, 0.70423929, 0.72666444, 0.74908959, 0.77151474, 0.79393989, 0.81636504, 0.83879019, 0.86121534, 0.88364049, 0.90606564, 0.92849079, 0.95091593, 0.97334108, 0.99576623, 1.01819138, 1.04061653, 1.06304168, 1.08546683, 1.10789198, 1.13031713, 1.15274228, 1.17516743, 1.19759257, 1.22001772, 1.24244287, 1.26486802, 1.28729317, 1.30971832, 1.33214347, 1.35456862, 1.37699377, 1.39941892, 1.42184407, 1.44426921, 1.46669436, 1.48911951, 1.51154466, 1.53396981, 1.55639496, 1.57882011, 1.60124526, 1.62367041, 1.64609556, 1.66852071, 1.69094585, 1.713371 , 1.73579615, 1.7582213 , 1.78064645, 1.8030716 , 1.82549675, 1.8479219 , 1.87034705, 1.8927722 , 1.91519735, 1.93762249, 1.96004764, 1.98247279, 2.00489794]), <a list of 100 Patch objects>)