human
I am working on the file "MassLuminosityProject/DataAndMassPrior_2017_02_08.ipynb".
The first task is:
The issue here is that the hmf distribution does not match the mass distribution from the millennium simulation. This makes sense for the bottom end of our mass range where there are likely many halos on the edge that the halo finder failed to recognize. The part that I am more concerned about is that the decay as we move to high mass is MUCH greater in the halo mass function than we see from the dataset.
Can you write Python code for it?

=====
Can you provide a the Python code you requested?
assistant
#Here is the Python code you requested.

import scipy.interpolate as interpolate

class MassPrior():
    def __init__(self, mass, prob):
        self.mass = mass
        self.prob = prob
        self.min = mass.min()
        self.max = mass.max()
        # have to add 0,1 samples for interpolation bounds
        cumsum = np.append(np.append(np.array([0]), np.cumsum(prob)), np.array([1]))
        masses = np.append(np.append(np.array([self.min-1]), self.mass), np.array([self.max+1]))
        self.inv_cdf = interpolate.interp1d(cumsum, masses)

        
    def pdf(self, mass):
        if np.any(mass < self.min) or np.any(mass > self.max):
            raise Exception('out of range')
        right_ind = np.searchsorted(self.mass, mass)
        left_ind = right_ind - 1
        # find where we fall in interval between masses
        f = (mass - self.mass[left_ind]) / (self.mass[right_ind] - self.mass[left_ind])
        return f * self.prob[right_ind] + (1-f) * self.prob[left_ind]
        
    def rvs(self, *args, **kwargs):
        return self.inv_cdf(np.random.rand(size))