Import NumPy and SciPy (not needed when using --pylab)
In [ ]:
%pylab inline
import scipy as sp
Load data from file
In [ ]:
zz = np.loadtxt('wiggleZ_DR1_z.dat',dtype='float'); # Load WiggleZ redshifts
Check bounds
In [ ]:
np.min() # Check bounds
In [ ]:
np.max()
Construct histogram from data
There are several histrogram commands: hist() will be fine here, but note the syntax below. Also note that the bin edges are returned, so that there will nbins+1 of these.
In [ ]:
nbins = 50; # Is this a good choice?
In [ ]:
n, bins, patches = hist() # With hist, one needs to (spuriously) request the patch objects as well
In [ ]:
x = bins[0:nbins] + (bins[2]-bins[1])/2; # Convert bin edges to centres, chopping the last
Interpolate histogram output -> p(z); n.b. that you can also use numerical quadrature to get $P(z)$ directly.
In [ ]:
# Import the function you need
from scipy.interpolate import interp1d
In [ ]:
# Build an interpolation function for p(z) that accepts an arbitrary redshift z
In [ ]:
z = linspace(0,2,100); plot(z,p(z)) # Test your interpolation function out
Use numerical integration to get $P(z) = \int_0^\infty p(z') dz'$
In [ ]:
# Import the function you need
from scipy import integrate
In [ ]:
Pz = lambda : ... # Use integrate inside a lambda function to define P(z)?
In [ ]:
total = Pz(5) # Get normalisation constant by evaluating P(z->\infty)
In [ ]:
total # Check that this worked