mrpy¶
A Python package for calculations with the MRP parameterisation of the Halo Mass Function.
See Murray, Robotham, Power (2018) for more details on what the MRP is.
Installation¶
>> pip install mrpy.
This should install the required dependencies automatically.
Note, to use the MCMC fitting features, emcee is needed. This is not installed automatically.
To get the bleeding edge, use pip install git+git://github.com/steven-murray/mrpy.git.
Getting Started¶
There’s a lot of things that you can do with mrpy. What you require will depend on the problem at hand. We recommend looking at some of the examples, and the API itself for how to use the code.
Features¶
With mrpy you can:
- Calculate basic statistics of the truncated generalised gamma distribution (TGGD) with the TGGD class: mean, mode, variance, skewness, pdf, cdf, generate random variates etc.
- Generate MRP quantities with the MRP class: differential number counts, cumulative number counts, various methods for generating normalisations.
- Generate the MRP-based halo mass function as a function of physical parameters via the mrp_b13 function.
- Fit MRP parameters to data in the form of arbitrary curves with the get_fit_curve function.
- Fit MRP parameters to data in the form of a sample of variates with the SimFit class: simulation data is supported with extra efficiency, simulation suites fitted simultaneously is also supported, arbitrary priors on parameters, log-normal uncertainties on variates supported.
- Calculate analytic hessians, jacobians at any point (including the solution of a fit).
- Use alternate parameterisations of the same form via the reparameterise module.
- Work with a special entirely analytic model to understand the effects of various parameters in the analytic_model module.
Examples¶
There are several examples featured in the docs/examples directory of the github repository. These can also be found
in the official documentation.
Acknowledging¶
If you use this code in your work, please cite Murray, Robotham, Power (2018) and/or http://ascl.net/1802.015 (whichever is more appropriate). Also consider starring/following the repo on github so we know how much it is being used. We would also love any input to the code!
Contents¶
Examples¶
To help get you started using mrpy, we’ve compiled a few examples. In fact, many of these examples are used explicitly
to produce the figures which appear in the MRP paper (Murray, Robotham, Power 2018). Other simple examples can be found in
the API documentation for each object.
So, what would you like to learn?
How to get started¶
This series of examples shows the very basics of how to get started with MRPY, using the different functionality. These examples aren’t “real world” ones, just toy ones to show the basic idea.
In [1]:
# General Imports
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
Core Functionality¶
Core functionality (i.e. calculation of the MRP function given input
parameters, plus some other functions useful for normalising) is in the
core module:
In [2]:
from mrpy import dndm # resides in the base.core module, but imported into top-level namespace
m = np.logspace(10,15,500) # Create an array of masses
dn = dndm(m, logHs=14.0, alpha=-1.8, beta=0.85) # Create the MRP mass function
plt.plot(m,dn,lw=2)
plt.xscale('log')
plt.yscale('log')
Pure Stats¶
If you don’t care so much about the fact that the MRP is good for halo
mass functions (or don’t know what a halo mass function is…), but want
to use the statistical distribution, you’ll want the stats module.
It contains an object called TGGD (short for Truncated Generalised
Gamma Distribution), which has many statistical quantities and methods
available (such as producing random variates, mean, mode etc.)
In [3]:
from mrpy import TGGD
tggd = TGGD(scale=1e14,a=1.0,b=0.85,xmin=1e10)
print "Mean: ", tggd.mean # Mean of the distribution
print "Mode: ", tggd.mode # Mode of the distirbution
print "Variance: ", tggd.variance # Variance of the distribution
print "Mean of sample: ", np.mean(tggd.rvs(1e5)) #Produce 1e5 random variates and take the mean
plt.plot(m, tggd.pdf(m)) #Plot the PDF
plt.xscale('log')
Mean: 2.84867015869e+14
Mode: 1.21070004341e+14
Variance: 4.79705281452e+28
Mean of sample: 2.83287193279e+14
Physical Dependence¶
The physical_dependence module contains a counterpart to the basic
dndm function, called mrp_b13, which returns the best-fit MRP
according to input physical variables (redshift, matter density, rms
mass variance). These are derived from fits to the theoretical mass
function of Behroozi+2013.
In [4]:
from mrpy import mrp_b13
dndm_z1 = mrp_b13(m,z=1) # HMF at redshift 1
dndm_z0 = mrp_b13(m,sigma_8=0.85) # HMF at redshift 0 but sigma_8=0.85
plt.plot(m,dndm_z1,label="z=1")
plt.plot(m,dndm_z0,label="z=0")
plt.legend(loc=0)
plt.xscale('log')
plt.yscale('log')
Fitting MRP¶
The fit_curve module contains routines to fit the MRP to
binned/curve data. This can be a theoretical curve, or binned halos (or
other variates). There are several options available, and the gradient
of the objective function is specified analytically to improve
performance. See Murray, Robotham, Power, 2016 (in prep.) for more
details.
In [5]:
from mrpy.fitting.fit_curve import get_fit_curve
dn = dndm(m,logHs=14.0,alpha=-1.9,beta=0.75,norm=1)
result, curve_obj = get_fit_curve(m,dn,[14.5, -1.8, 0.85,0.],
bounds = [[None,None], [-2,-1.5], [0.2,1.5], [None,None]]) #The bounds argument is very important at this stage.
print result
plt.plot(curve_obj.logm,curve_obj.dndm()/dn-1,lw=2)
fun: 1.6704623362075916e-12
hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
jac: array([ -8.60315661e-05, -1.13296916e-04, 6.03274702e-05,
-3.79420248e-06])
message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
nfev: 28
nit: 25
status: 0
success: True
x: array([ 1.40000001e+01, -1.90000003e+00, 7.50000058e-01,
-5.62234817e-07])
Out[5]:
[<matplotlib.lines.Line2D at 0x7f072e383dd0>]
In the previous example, we simply fit the four MRP parameters to the input curve. Options can be specified to constrain the normalisation via the known mean density of the Universe (see MRP, 2017 Appendix C.1 for details).
To fit actual samples of halos, use the fit_sample module. At this
stage, only samples with no measurement uncertainties, and a constant
volume (per subsample) are supported. Usefully, this covers the case of
output halos from simulations. Within this context, either a
downhill-gradient method or MCMC can be used.
We hope to provide more general fitting scenarios via configuration with other packages in the future.
An example of using the downhill method is as follows:
In [6]:
from mrpy.base.stats import TGGD
from mrpy.fitting.fit_sample import SimFit
# Create some mock data to fit
r = TGGD(scale=1e14,a=-1.8,b=1.0,xmin=1e12).rvs(1e5)
r = np.sort(r)
# Create fit object, specifying parameter bounds
fitobj = SimFit(r,hs_bounds=(12,16),alpha_bounds=(-1.99,-1.6),beta_bounds=(0.5,1.5))
# Run downhill gradient method
res, obj = fitobj.run_downhill()
# Print the resulting parameters
print "Resulting Parameters: ", res.x
# The "obj" returned is a list of PerObjLike objects defined at the result, containing lots of methods and quantities:
print "Hessian at solution: ", obj[0].hessian
print "Covariance at solution: ", obj[0].cov
# Plot the mass function from obj (the logm contains all the masses from the fit)
plt.plot(obj[0].logm,obj[0].dndm(log=True),lw=2)
Resulting Parameters: [ 13.9900839 -1.79700413 0.94765685 -24.43903982]
Hessian at solution: [[-1838656.7724181 1485394.57073756 -508029.87769594 -427802.81517913]
[ 1485394.57073756 -1328344.72287371 412360.36967558 351970.2695528 ]
[ -508029.87769594 412360.36967558 -142628.15625157 -118826.3858363 ]
[ -427802.81517913 351970.2695528 -118826.3858363 -100000.08553355]]
Covariance at solution: [[ 0.00077558 -0.00025638 0.00122131 -0.00567158]
[-0.00025638 0.00010004 -0.00046645 0.00200319]
[ 0.00122131 -0.00046645 0.00287933 -0.01028797]
[-0.00567158 0.00200319 -0.01028797 0.04354857]]
Out[6]:
[<matplotlib.lines.Line2D at 0x7f072cbebd50>]
Fit MRP parameters to model data¶
In this example, we do something very simple – fit the MRP parameters using a theoretically produced HMF. This might be one of the first things you’d want to do with the MRP. In addition to the simple fit, we’ll also change the truncation scale, to see how the fit performs over different mass ranges. Furthermore, we’ll change the redshift and halo overdensity to make sure the fit performs well in all cases.
The resulting plots appear as figures 1,2 in Murray, Robotham, Power (2018)
In [28]:
#Import necessary libraries
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import MaxNLocator
from mrpy.fitting.fit_curve import get_fit_curve
from os.path import join, expanduser
#We'll use the hmf package to produce the theoretical HMF
from hmf import MassFunction
from hmf.cosmo import Planck13
In [17]:
fig_folder = expanduser("~")
First, we produce a hmf model consistent with the Planck 2013 data, and resolved enough to produce a high-quality fit:
In [29]:
hmf = MassFunction(hmf_model="Tinker08",Mmin=5,Mmax=18,lnk_min=-13,lnk_max=13,
transfer_model="CAMB",
dlnk=0.01,
sigma_8=0.829,n=0.9603,dlog10m=0.1,cosmo_model=Planck13)
Here’s the important part: actually fitting the data.
Fit across redshift and overdensity¶
In [13]:
# Create the lists that we'll use to save the results
res = [0]*64
obj = [0]*64
max_dev = np.zeros(64)
rms_dev = np.zeros(64)
dndms_whole = [0]*64
ms_whole = [0]*64
# 4 different truncation masses
sigmins = [4,3,2,1.5]
mmins = [1e10,1e11,1e12,1e13]
deltahs = [200.0,400.0,800.0,1600.0]
zs = [0.0,0.5,1.0,2.0]
dndms = [0]*64
ii = 0
for i,z in enumerate(zs):
for j,deltah in enumerate(deltahs):
hmf.update(z=z,delta_h=deltah)
for k,sigmin in enumerate(sigmins):
# Get theoretical data
mask = np.logical_and(hmf.sigma < sigmin, hmf.sigma > 0.5)
dn = hmf.dndm[mask]
m = hmf.m[mask]
# Fit the MRP in the *simplest* way possible.
res[ii], obj[ii] = get_fit_curve(m, dn,
x0 = [14.4,-1.9,0.8,-43],
bounds=[[None,None],[-2.5,-1.5],[0.2,None],[None,None]],
jac=False)
max_dev[ii] = np.abs(obj[ii].dndm()/dn-1).max()
rms_dev[ii] = np.sqrt(np.mean((obj[ii].dndm()/dn-1)**2))
dndms[ii] = dn
dndms_whole[ii] = hmf.dndm
ms_whole[ii] = hmf.m
#print z, deltah, "%1.0e : "%mmin, res[ii].x, max_dev[ii], rms_dev[ii] # the actual result values
ii += 1
Next the boring stuff… setting up and plotting the figure.
In [18]:
# Create the figure object
xmin_plot = 1e7
xmax_plot = 3e15
ymin_plot = -0.08
ymax_plot = 0.08
fig,ax = plt.subplots(4,4,sharex=True,sharey=True,squeeze=True,figsize=(12,8),
subplot_kw={"xscale":"log","ylim":(ymin_plot,ymax_plot),
"xlim":(xmin_plot, xmax_plot)})
# Contract the space a bit
plt.subplots_adjust(wspace=0.08,hspace=0.10)
# Plot the fitted data.
# Note that 'obj' contains lots of quantities of interest, not least of which is a method
# to calculate dn/dm!
ii = 0
for i,z in enumerate(zs):
for j,deltah in enumerate(deltahs):
ax[i,j].text(2*xmin_plot,0.053, "RMS(%)" + r"$\sim %1.1f - %1.1f$"%(rms_dev[ii:ii+4].min()*100,rms_dev[ii:ii+4].max()*100),fontsize=13)
ax[i,j].text(2*xmin_plot,-0.06,r"$z=%s,\ \Delta_h=%s$"%(z,deltah),fontsize=13)
for k,mmin in enumerate(mmins):
# Background grey scaled MF
if k==0:
mask = np.logical_and(ms_whole[ii]>xmin_plot, ms_whole[ii]<xmax_plot)
dndm = np.log10(dndms_whole[ii][mask])
if ii==0:
dndm_range = np.max(dndm) - np.min(dndm)
dndm *= (ymax_plot - ymin_plot)/dndm_range
dndm += ymax_plot - np.max(dndm)
ax[i,j].plot(ms_whole[ii][mask], dndm, color='grey', alpha=0.4,lw=3)
# Plot each iteration
ax[i,j].plot(obj[ii].m,obj[ii].dndm()/dndms[ii]-1,lw=2)
# Modify the ticks for prettiness
ax[i,j].tick_params(axis='both', which='major', labelsize=11)
ax[i,j].tick_params(axis='both', which='major', labelsize=11)
ax[i,j].yaxis.set_major_locator(MaxNLocator(6))
ii += 1
fig.text(0.5, 0.04, r"Mass, $h^{-1}M_\odot$",fontsize=15, ha='center', va='center')
fig.text(0.06, 0.5, 'Relative Residuals', fontsize=15,ha='center', va='center', rotation='vertical')
#Save image
if fig_folder:
fig.savefig(join(fig_folder,"comparison_tinker.pdf"))
Fit against different cosmology and halo type¶
In [22]:
hmf = MassFunction(hmf_model="Warren",Mmin=5,Mmax=18,lnk_min=-13,lnk_max=13,
transfer_model="CAMB",
dlnk=0.01,
sigma_8=0.829,n=0.9603,dlog10m=0.1,cosmo_model=Planck13)
In [23]:
# Create the lists that we'll use to save the results
res = [0]*64
obj = [0]*64
max_dev = np.zeros(64)
rms_dev = np.zeros(64)
# 4 different truncation masses
sigmins = [4,3,2,1.5]
oms = [0.2, 0.25, 0.3, 0.35]
s8s = [0.7, 0.8, 0.9, 1.0]
dndms = [0]*64
ii = 0
for i,s8 in enumerate(s8s):
for j,om in enumerate(oms):
hmf.update(cosmo_params={"Om0":om}, sigma_8 = s8)
for k,sigmin in enumerate(sigmins):
# Get theoretical data
mask = np.logical_and(hmf.sigma < sigmin, hmf.sigma > 0.5)
dn = hmf.dndm[mask]
m = hmf.m[mask]
# Fit the MRP in the *simplest* way possible.
res[ii], obj[ii] = get_fit_curve(m, dn,
x0 = [14.4,-1.9,0.8,-43],
bounds=[[None,None],[-2.5,-1.5],[0.2,None],[None,None]],
jac=False)
max_dev[ii] = np.abs(obj[ii].dndm()/dn-1).max()
rms_dev[ii] = np.sqrt(np.mean((obj[ii].dndm()/dn-1)**2))
dndms[ii] = dn
dndms_whole[ii] = hmf.dndm
ms_whole[ii] = hmf.m
#print z, deltah, "%1.0e : "%mmin, res[ii].x, max_dev[ii], rms_dev[ii] # the actual result values
ii += 1
In [25]:
xmin_plot = 1e7
xmax_plot = 3e15
ymin_plot = -0.08
ymax_plot = 0.08
fig,ax = plt.subplots(4,4,sharex=True,sharey=True,squeeze=True,figsize=(12,8),
subplot_kw={"xscale":"log","ylim":(ymin_plot,ymax_plot),
"xlim":(xmin_plot, xmax_plot)})
# Contract the space a bit
plt.subplots_adjust(wspace=0.08,hspace=0.10)
# Plot the fitted data.
# Note that 'obj' contains lots of quantities of interest, not least of which is a method
# to calculate dn/dm!
ii = 0
for i,s8 in enumerate(s8s):
for j,om in enumerate(oms):
ax[i,j].text(2*xmin_plot,0.053, "RMS(%)" + r"$\sim %1.1f - %1.1f$"%(rms_dev[ii:ii+4].min()*100,rms_dev[ii:ii+4].max()*100),fontsize=13)
ax[i,j].text(2*xmin_plot,-0.06,r"$\Omega_m=%s,\ \sigma_8=%s$"%(om,s8),fontsize=13)
for k,mmin in enumerate(mmins):
if k==0:
mask = np.logical_and(ms_whole[ii]>xmin_plot, ms_whole[ii]<xmax_plot)
dndm = np.log10(dndms_whole[ii][mask])
if ii==0:
dndm_range = np.max(dndm) - np.min(dndm)
dndm *= (ymax_plot - ymin_plot)/dndm_range
dndm += ymax_plot - np.max(dndm)
ax[i,j].plot(ms_whole[ii][mask], dndm, color='grey', alpha=0.4,lw=3)
# Plot each iteration
ax[i,j].plot(obj[ii].m,obj[ii].dndm()/dndms[ii]-1,lw=2)
# Modify the ticks for prettiness
ax[i,j].tick_params(axis='both', which='major', labelsize=11)
ax[i,j].tick_params(axis='both', which='major', labelsize=11)
ax[i,j].yaxis.set_major_locator(MaxNLocator(6))
ii += 1
fig.text(0.5, 0.04, r"Mass, $h^{-1}M_\odot$",fontsize=15, ha='center', va='center')
fig.text(0.06, 0.5, 'Relative Residuals', fontsize=15,ha='center', va='center', rotation='vertical')
#Save image
if fig_folder:
fig.savefig(join(fig_folder,"comparison_warren.pdf"))
Fit MRP parameters to a suite of simulation data simultaneously¶
In this example, we grab haloes from the publicly available \(\nu^2\)GC simulation suite and show how MRP can be fit to the haloes of 4 simulations simultaneously. In this case, the 4 simulations have different box sizes, so they probe different parts of the mass function more or less effectively. By combining them, we can get a good handle on a wide range of the mass function.
Do note that this example is not quick. It takes a while to get the data, let alone run the MCMC on it. You may want to generate some smaller fake datasets to have a play.
The plots from this example are used in MRP as Figures 3 and 4.
In [1]:
# General imports
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pickle
from os.path import join, splitext, exists, expanduser
# Mrpy imports
from mrpy.fitting.fit_sample import SimFit
from mrpy import MRP
from chainconsumer import ChainConsumer
In [2]:
fig_folder = expanduser("~")
data_folder = expanduser("~")
Preparing the data¶
First you’ll need to get the data. These are the files you’ll need to download (beware, at least one of them is 12Gb alone):
http://www2.ccs.tsukuba.ac.jp/Astro/Members/ishiyama/nngc/Data/n2gc-m_z0.fof.bz2
http://www2.ccs.tsukuba.ac.jp/Astro/Members/ishiyama/nngc/Data/n2gc-h1_z0.fof.bz2
http://www2.ccs.tsukuba.ac.jp/Astro/Members/ishiyama/nngc/Data/n2gc-m_z0.rockstar.bz2
Then unzip them. NOTE: you don’t need to run this section if you’ve already got the data and compactified it
First of all, we need to pare down the huge files. We can do this in a few ways:
- We only care about the Mass column, so we can delete everything else
- We only care about unique halo masses (and the quantity of each), so we can “tabulate” the data
- We keep only haloes with 40 or more particles (this limit is taken from the paper accompanying the catalogue)
These operations reduce the file size by about a factor of 100-1000, and make the subsequent MCMC runs much faster.
Something else to consider is that the fastest way to read in the data files and reduce them is to do it in one big chunk with numpy. However, this takes a lot of memory. So instead we read them line by line.
With these considerations, we implement the following functions.
In [3]:
def strip_and_compress(fname,fout,mpart=None,Nmin=0,Nmax=np.inf, force=False):
unique_masses = {}
ftype = splitext(fname)[1][1:]
if not force and exists(fout):
return
with open(fname) as fin:
for line in fin:
l = line.strip()
# Skip comments
if l.startswith("#"):
continue
else:
if ftype=="fof":
npart = int(l.split()[-1])
elif ftype=="rockstar":
npart = int(l.split()[7])
# Reject the entry if it is less than Nmin
if npart < Nmin:
continue
elif npart > Nmax:
continue
# Calculate the mass of the halo
if ftype=="fof":
mvir = mpart * npart
elif ftype=="rockstar":
mvir = float(l.split()[21]) # Corresponds to M200b
# Add it to the final unique mass dict
if mvir in unique_masses:
unique_masses[mvir] += 1
else:
unique_masses[mvir] = 1
# Convert the dict of values into a 2D array of masses and number of occurrences
out = np.array([[k,v] for k,v in unique_masses.iteritems()])
print "Compressed {} to {} percent".format(fname,100*len(out[:,1])/sum(out[:,1]))
# Save the data to a table file
np.savetxt(fout,out)
Now actually do the stripping and compressing of the files. We save the data in new files with an appended “.compact”. Note also we limit the size of the halos, to be in line with the quoted values from the I15 paper. There is in fact at least 1 outlier beyond these limits.
In [76]:
FORCE = False
strip_and_compress(join(data_folder,"n2gc-h1_z0.fof"),
join(data_folder, "n2gc-h1_z0.fof.compact"),2.75e7,100,17476256)
strip_and_compress(join(data_folder,"n2gc-m_z0.fof"),
join(data_folder, "n2gc-m_z0.fof.compact"),2.2e8,100,12120576)
strip_and_compress(join(data_folder,"n2gc-m_z0.rockstar"),
join(data_folder, "n2gc-m_z0.rockstar.compact"),force=FORCE,
Nmin=100)
strip_and_compress(join(data_folder,"n2gc-h1_z0.rockstar"),
join(data_folder, "n2gc-h1_z0.rockstar.compact"), force=FORCE,
Nmin=100)
Compressed /home/steven/Documents/DataSets/n2gc/n2gc-m_z0.rockstar to 0.219503952011 percent
Compressed /home/steven/Documents/DataSets/n2gc/n2gc-h1_z0.rockstar to 0.949670196418 percent
First up, read in the compact data we just created.
In [4]:
# Read in the data from file
def get_raw_data(folder, sims=['h1','m'],ftype="fof",mmin=None,mmax=None):
m = []
nm = []
for sim in sims:
data = np.genfromtxt(join(folder,"n2gc-{}_z0.{}.compact".format(sim,ftype)))
m.append(data[:,0])
nm.append(data[:,1])
if mmin is not None:
for i,mm in enumerate(mmin):
nm[i] = nm[i][m[i]>mm]
m[i] = m[i][m[i]>mm]
if mmax is not None:
for i,mm in enumerate(mmax):
nm[i] = nm[i][m[i]<mm]
m[i] = m[i][m[i]<mm]
return m,nm
We read in both FOF and SO halos with similar parameters, and store
everything in the data dictionary.
In [77]:
data = {'fof':{}, "so":{}}
# FOF halos
data['fof']['m'], data['fof']['nm'] = get_raw_data(data_folder, ['h1','m'],
mmin=[2.75e9,2.2e10],mmax=[2e13,7e14])
data['fof']['weights'] = [data['fof']['nm'][0]/140.0**3, data['fof']['nm'][1]/560.0**3]
print "Total number of FOF haloes: ", np.sum([np.sum(x) for x in data['fof']['nm']])
print "Total number of *unique* FOF haloes: ", np.sum([len(x) for x in data['fof']['m']])
print "-"*40
# SO halos
data['so']['m'], data['so']['nm'] = get_raw_data(data_folder, ['h1','m'], ftype='rockstar',
mmin=[2.75e9,2.2e10],mmax=[2e13,7e14])
data['so']['weights'] = [data['so']['nm'][0]/140.0**3, data['so']['nm'][1]/560.0**3]
print "Total number of SO haloes: ", np.sum([np.sum(x) for x in data['so']['nm']])
print "Total number of *unique* SO haloes: ", np.sum([len(x) for x in data['so']['m']])
Total number of FOF haloes: 24640920.0
Total number of *unique* FOF haloes: 141278
----------------------------------------
Total number of SO haloes: 24245670.0
Total number of *unique* SO haloes: 109419
Running the fits¶
We’ll run the fits with the emcee package (via a routine built in to
mrpy), but also with an optimization solver. The in-built function
is able to utilise the tabulation of data we have performed already, and
can do the suites simultaneously.
In [78]:
# Create the fitting class instance. This will have uniform priors.
fitobj_fof = SimFit(data['fof']['m'],data['fof']['nm'],
V=[140.0**3,560.0**3],
alpha_bounds = (-1.99,-1.5), hs_bounds=(12,16),
beta_bounds=(0.2,1.5),lnA_bounds=(-50,-10))
fitobj_so = SimFit(data['so']['m'],data['so']['nm'],
V=[140.0**3,560.0**3],
alpha_bounds = (-1.99,-1.5), hs_bounds=(12,16),
beta_bounds=(0.2,1.5),lnA_bounds=(-50,-10))
In [79]:
# We don't use these, but they can be useful if something goes wrong.
downhill_res_fof = fitobj_fof.run_downhill(lnA0=-40.0)
downhill_res_so = fitobj_so.run_downhill(lnA0=-40.0)
In [80]:
# Run the mcmc.
# We set 300 chains to warmup, but we can extend this later if we need to manually.
# Also, we start the chains in a small ball around the best (downhill) optimization solution using opt_init=True.
#fitobj_fof.run_mcmc(nchains=50,warmup=200,iterations=500,opt_init=True,threads=8)
fitobj_so.run_mcmc(nchains=50,warmup=200,iterations=500,opt_init=True,threads=8)
Out[80]:
<emcee.ensemble.EnsembleSampler at 0x7fc6704bd350>
First off we want to look at a few key diagnostics of the chains to check whether everything’s okay.
In [81]:
print "Acceptance fraction for FOF (min, max, mean): ", fitobj_fof.mcmc_res.acceptance_fraction.min(), fitobj_fof.mcmc_res.acceptance_fraction.max(), fitobj_fof.mcmc_res.acceptance_fraction.mean()
print "Acceptance fraction for SO (min, max, mean): ", fitobj_so.mcmc_res.acceptance_fraction.min(), fitobj_so.mcmc_res.acceptance_fraction.max(), fitobj_so.mcmc_res.acceptance_fraction.mean()
Acceptance fraction for FOF (min, max, mean): 0.518 0.622 0.57288
Acceptance fraction for SO (min, max, mean): 0.542 0.648 0.58912
These acceptance fractions are somewhat high, but probably okay. We’ll check burnin as well soon.
In [5]:
def gelman_rubin(chain):
ssq = np.var(chain, axis=1, ddof=1)
W = np.mean(ssq, axis=0)
thb = np.mean(chain, axis=1)
thbb = np.mean(thb, axis=0)
m = chain.shape[0]
n = chain.shape[1]
B = n / (m - 1) * np.sum((thbb - thb)**2, axis=0)
var_th = (n - 1.) / n * W + 1. / n * B
R = np.sqrt(var_th / W)
return R
In [82]:
ChainConsumer().add_chain(fitobj_fof.mcmc_res.chain.reshape((-1,4)), walkers = 50).diagnostic.gelman_rubin(threshold=0.1)
ChainConsumer().add_chain(fitobj_so.mcmc_res.chain.reshape((-1,4)), walkers = 50).diagnostic.gelman_rubin(threshold=0.1)
Gelman-Rubin Statistic values for chain 0
Param 0: 1.04796 (Passed)
Param 1: 1.04777 (Passed)
Param 2: 1.09876 (Passed)
Param 3: 1.05742 (Passed)
Gelman-Rubin Statistic values for chain 0
Param 0: 1.06547 (Passed)
Param 1: 1.06030 (Passed)
Param 2: 1.06584 (Passed)
Param 3: 1.06507 (Passed)
Out[82]:
True
We see that the chains have converged (R < 1.1).
Since the fitting takes some time, we save the main results, i.e. the chain, to file here so that we can begin again at any time without running the MCMC. Thus the following analysis only uses the chains as written to file, rather than the full fit objects just created.
In [83]:
np.savez("n2gc_analysis/n2gc_mcmc_chain_fof",chain=fitobj_fof.mcmc_res.chain)
np.savez("n2gc_analysis/n2gc_mcmc_chain_so",chain=fitobj_so.mcmc_res.chain)
In [6]:
chain_so = np.load("n2gc_analysis/n2gc_mcmc_chain_so.npz")['chain']
chain_fof = np.load("n2gc_analysis/n2gc_mcmc_chain_fof.npz")['chain']
Analysis¶
The first thing we might want to do with each fit is to check its traceplot, and determine if the burnin was sufficient.
In [7]:
def traceplot(keys,chains):
f, ax = plt.subplots(len(keys), 1, sharex=True, figsize=(8, 2.5 * len(keys)))
for i, (key, chain) in enumerate(zip(keys,chains.T)):
ax[i].plot(chain, color="black", alpha=0.2)
ax[i].set_ylabel(key,fontsize=16)
f.subplots_adjust(hspace=0)
plt.setp([a.get_xticklabels() for a in f.axes[:-1]], visible=False)
return f
Plot both cases:
In [44]:
fig0 = traceplot([r"$ \log \mathcal{H}_\star$",r"$\alpha$",r'$\beta$',r"$\ln A$"],chain_fof)
plt.show()
It seems like a reasonable burn-in time has been met for the FOF halos, so we’re happy we trust our sample.
In [86]:
fig0 = traceplot([r"$ \log \mathcal{H}_\star$",r"$\alpha$",r'$\beta$',r"$\ln A$"],chain_so)
plt.show()
The SO halos however have moved significantly, possible due to a poor downhill-gradient optimization. We remove the firs parts of the chain:
In [87]:
chain_so = chain_so[:,100:,:]
We’d like to know the pure basic results: mean, median, mode, standard deviation etc.
In [51]:
print "Mean: ", np.mean(chain_fof,axis=(0,1))
print "Median: ", np.median(chain_fof,axis=(0,1))
print "Mode: ", fitobj_fof.mcmc_res.flatchain[np.argmax(fitobj_fof.mcmc_res.flatlnprobability),:]
print "Std Dev.: ", np.std(chain_fof,axis=(0,1))
print "Covariance: ", np.cov(chain_fof.reshape((-1,4)).T)
print "Relative Uncertainty: ", np.std(chain_fof,axis=(0,1))*100/np.mean(chain_fof,axis=(0,1))
from mrpy.base.core import log_mass_mode
print "Log Mass Mode: ", np.log10(log_mass_mode(*np.mean(chain_fof[:,:,:3],axis=(0,1))))
Mean: [ 14.51918594 -1.89972215 1.19198794 -44.41676165]
Median: [ 14.51921005 -1.89971928 1.19135641 -44.41679216]
Mode: [ 14.51964269 -1.89972233 1.19189829 -44.41885557]
Std Dev.: [ 0.00715219 0.00016809 0.0194396 0.03367548]
Covariance: [[ 5.11558477e-05 -1.68702393e-07 -1.92940482e-05 -2.09231989e-04]
[ -1.68702393e-07 2.82549206e-08 -1.01344465e-06 1.83532519e-06]
[ -1.92940482e-05 -1.01344465e-06 3.77913140e-04 -2.41123515e-04]
[ -2.09231989e-04 1.83532519e-06 -2.41123515e-04 1.13408347e-03]]
Relative Uncertainty: [ 0.04926026 -0.00884806 1.63085537 -0.07581706]
Log Mass Mode: 13.617275083
In [88]:
print "Mean: ", np.mean(chain_so,axis=(0,1))
print "Median: ", np.median(chain_so,axis=(0,1))
print "Mode: ", fitobj_so.mcmc_res.flatchain[np.argmax(fitobj_so.mcmc_res.flatlnprobability),:]
print "Std Dev.: ", np.std(chain_so,axis=(0,1))
print "Covariance: ", np.cov(chain_so.reshape((-1,4)).T)
print "Relative Uncertainty: ", np.std(chain_so,axis=(0,1))*100/np.mean(chain_so,axis=(0,1))
from mrpy.base.core import log_mass_mode
print "Log Mass Mode: ", np.log10(log_mass_mode(*np.mean(chain_so[:,:,:3],axis=(0,1))))
Mean: [ 14.10470648 -1.88130297 0.74993505 -42.00895938]
Median: [ 14.10465981 -1.88130181 0.7497424 -42.00886146]
Mode: [ 14.10427937 -1.88133126 0.75064434 -42.00835854]
Std Dev.: [ 0.00592467 0.00024254 0.00635023 0.02659803]
Covariance: [[ 3.51034875e-05 -4.02509937e-07 -6.81369166e-06 -1.46285387e-04]
[ -4.02509937e-07 5.88305367e-08 -9.58809971e-07 3.51176102e-06]
[ -6.81369166e-06 -9.58809971e-07 4.03274777e-05 -3.28788684e-05]
[ -1.46285387e-04 3.51176102e-06 -3.28788684e-05 7.07490346e-04]]
Relative Uncertainty: [ 0.04200493 -0.01289234 0.84677108 -0.06331513]
Log Mass Mode: 13.0371689695
This produces a “corner” plot which shows the covariance between parameters.
In [53]:
c = ChainConsumer().add_chain(chain_fof.reshape((-1,4)),
parameters=[r'$h_\star$',r'$\alpha$',r"$\beta$",r'$\ln A$'],
walkers=50)
fig = c.plot(figsize="PAGE")
if fig_folder:
fig.savefig(join(fig_folder,"n2gc_triangle.pdf"))
WARNING:chainconsumer.chain:This method is deprecated. Please use chainConsumer.plotter.plot instead
In [89]:
c = ChainConsumer().add_chain(chain_so.reshape((-1,4)),
parameters=[r'$h_\star$',r'$\alpha$',r"$\beta$",r'$\ln A$'],
walkers=50)
fig = c.plot(figsize="PAGE")
WARNING:chainconsumer.chain:This method is deprecated. Please use chainConsumer.plotter.plot instead
Importantly, we want to check if the actual results look good against the data, when binned.
In [90]:
# A function to create histograms from raw masses, and conver them to dn/dm.
# It also sets edge values in which a whole bin is not sampled to nan for visual purposes.
def bin_masses(masses,nm, V, bins=50):
hist, edges = np.histogram(np.log10(masses), bins,weights=nm)
centres = (edges[1:] + edges[:-1]) / 2
dx = centres[1] - centres[0]
dn = hist.astype("float") / (10 ** centres *float(V)* dx * np.log(10))#
poisson_error = np.sqrt(hist.astype("float"))/ (10 ** centres *float(V)* dx * np.log(10))#
try:
hist0 = np.where(hist != 0)[0][0]
dn[hist0] = 0
hist[hist0] = 0
poisson_error[hist0] = 0
except IndexError:
pass
try:
histN = np.where(hist != 0)[0][-1]
dn[histN] = 0
hist[histN] = 0
poisson_error[histN] = 0
except IndexError:
pass
dn[hist == 0] = np.nan
return centres, dn, hist, poisson_error
resids = {}
for jj, ftype in enumerate(['fof','so']):
resids[ftype] = {}
m,nm = data[ftype]['m'], data[ftype]['nm']
# Generate total density of each sim
resids[ftype]['rho'] = [np.sum(x*nx)/L**3 for x,nx,L in zip(m,nm,[140.0,560.0])]
# Calculate the total mmin and mmax for all sims in the suite
mmin = np.min([x.min() for x in m])
mmax = np.max([x.max() for x in m])
# Generate the bin structure
bins = np.linspace(np.log10(mmin), np.log10(mmax),50)
bin_centres = (bins[1:] + bins[:-1])/2
# Generate the dn/dm from the sims
resids[ftype]["dndm"] = []
resids[ftype]["hist"] = []
resids[ftype]["err"] = []
for mi,nmi,L in zip(m,nm,[140.0,560.0]):
_,dn,h_, err = bin_masses(mi,nmi,L**3,bins)
resids[ftype]["dndm"].append(dn)
resids[ftype]["hist"].append(h_)
resids[ftype]["err"].append(err)
# The final best-fit object.
parms = np.mean([chain_fof, chain_so][jj], axis=(0,1))
norm = parms[3] # downhill_res[0].x[3]
resids[ftype]['fit'] = MRP(logm = bin_centres,logHs=parms[0],alpha=parms[1],beta=parms[2],norm=norm)
Along with the best-fit MRP, we want to show the published mass function of the data, which we get from the hmf package.
In [57]:
from hmf import MassFunction
h = MassFunction(hmf_model="Ishiyama", cosmo_params={"Om0":0.31, "Ob0":0.048, "H0":68.0},
sigma_8=0.83, n=0.96,lnk_min=-15, lnk_max=15, dlnk=0.01,Mmin=bin_centres[0],Mmax=bin_centres[-1]+0.001,
dlog10m=bin_centres[1]-bin_centres[0])
Finally we draw the actual plot.
In [114]:
fig, ax = plt.subplots(1,2, figsize=(9,4), sharex=True, sharey=True,
subplot_kw={"xscale":'log', 'ylim':(-0.2,0.2)},
gridspec_kw={"wspace":0.05})
ftypes = ['fof','so']
for jj, ftype in enumerate(ftypes):
for i,(dn,hst,err, label,col) in enumerate(zip(resids[ftype]['dndm'],
resids[ftype]['hist'],
resids[ftype]['err'],
["H1","M"],
["C0",'C2'])):
fit = resids[ftype]['fit']
# Plot alternative type in grey
if jj==1:
dn_, err_ = resids[ftypes[(jj+1)%2]]['dndm'][i], resids[ftypes[(jj+1)%2]]['err'][i]
frac = (dn_/fit.dndm()) - 1
err_ = err_/fit.dndm()
mask = np.abs(frac)<0.3
#ax[jj].plot(10**bin_centres[mask],frac[mask],label=ftypes[(jj+1)%2].upper() if i else "",
# color='k',lw=2, alpha=0.1)
ax[jj].fill_between(10**bin_centres[mask],
frac[mask] - err_[mask], frac[mask]+err_[mask],
label=ftypes[(jj+1)%2].upper() if i else "",
color='k', alpha=0.2, facecolor=None, edgecolor=None, lw=0)
# Residuals to MRP. Mask trailing bits so that poisson noise doesn't dominate the view
frac = (dn/fit.dndm()) - 1
err_ = err/fit.dndm()
mask = np.abs(frac)<0.3
ax[jj].plot(10**bin_centres[mask],frac[mask],label=label if not jj else "",color=col,lw=2)
ax[jj].fill_between(10**bin_centres[mask],
frac[mask] - err_[mask], frac[mask]+err_[mask],
color=col, alpha=0.2)
if ftype=="fof":
frac = (dn/h.dndm) - 1
err_ = err/h.dndm
mask = np.abs(frac)<0.3
ax[jj].plot(10**bin_centres[mask],frac[mask],color=col,lw=2, ls='--')
ax[jj].fill_between(10**bin_centres[mask],
frac[mask] - err_[mask], frac[mask]+err_[mask],
color=col, alpha=0.2, hatch='/')
ax[jj].set_xlabel(r"Halo Mass, [$h^{-1}M_\odot$]",fontsize=15)
ax[jj].grid(True)
# Rsidual of Rockstar to MRP
#frac = dndm_rock/fit.dndm() -1
#plt.plot(10**bin_centres[np.abs(frac)<0.3],frac[np.abs(frac)<0.3], color="C3",label="SO")
# Legend item for I15 fit
ax[0].plot([0],[0],label="Residual to I15",ls="--",color="k")
ax[0].set_title("FOF Halos")
ax[1].set_title("SO Halos")
# PLOT STYLING
#ax[0].xscale('log')
#plt.grid(True)
#plt.ylim((-0.2,0.2))
#plt.ylim((-0.05,0.05))
ax[0].set_ylabel("Sim/Fit - 1",fontsize=15)
for jj in range(2):
ax[jj].legend(loc=0,ncol=2)
# Save for the paper!
if fig_folder:
plt.savefig(join(fig_folder,"n2gc_fof_simul.pdf"))
/home/steven/miniconda3/envs/mrpy/lib/python2.7/site-packages/ipykernel/__main__.py:31: RuntimeWarning: invalid value encountered in less
/home/steven/miniconda3/envs/mrpy/lib/python2.7/site-packages/ipykernel/__main__.py:41: RuntimeWarning: invalid value encountered in less
/home/steven/miniconda3/envs/mrpy/lib/python2.7/site-packages/ipykernel/__main__.py:19: RuntimeWarning: invalid value encountered in less
We notice that the residuals from MRP are very similar in magnitude to those from the full EPS-based fit, over a fairly wide range of masses. Note that it seems that the MRP will diverge more significantly below the mass threshold than the EPS fit. In any case, both diverge significantly less than the same simulation with haloes found with a spherical overdensity technique.
Explore an entirely analytic model which includes no Poisson error¶
In this example, we explore the usefulness of a completely analytic solution to an ideal “sample”. For details, see Appendix D.3.3 of the MRP paper. In brief, the “ideal sample” is composed of non-Poisson limited haloes extracted from a pure MRP distribution within some physical volume. As such, the solution is a priori the input parameters of the MRP distribution. Finding the covariance of those parameters is our task.
The whole framework of this problem has already been implemented in
mrpy.extra.analytic_model. We will use that framework to answer some
questions.
In [1]:
# Imports
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from mrpy.extra.analytic_model import IdealAnalytic
In [2]:
# Fiducial Parameters
beta = 0.75
alpha = -1.85
hs = 14.5
V=400.0**3 # Physical volume
# Some constants
log_mmin = np.linspace(11,15,10)
What is the expected variance of parameters versustruncation mass?¶
In [3]:
fig,ax = plt.subplots(3,1,sharex=True,sharey=False,figsize=(6,7),subplot_kw={"yscale":'log'},
gridspec_kw={"hspace":0.0})
K = IdealAnalytic(log_mmin=log_mmin,logHs=hs,alpha=alpha,beta=beta,lnA=0.0)
ax[0].plot(log_mmin-hs,np.sqrt(K.cov[:,0,0])/np.abs(hs),lw=2)
ax[1].plot(log_mmin-hs,np.sqrt(K.cov[:,1,1])/np.abs(alpha),lw=2)
ax[2].plot(log_mmin-hs,np.sqrt(K.cov[:,2,2])/np.abs(beta),lw=2)
#ax[0].legend(loc=0,ncol=2)
ax[0].set_ylabel(r"$\sigma_{\log \mathcal{H}_\star}/|\log \mathcal{H}_\star|$",fontsize=15)
ax[1].set_ylabel(r"$\sigma_\alpha/|\alpha|$",fontsize=15)
ax[2].set_ylabel(r"$\sigma_\beta/|\beta|$",fontsize=15)
ax[2].set_xlabel(r"$\log_{10}\left(m_{min}/\mathcal{H}_\star\right)$",fontsize=15)
Out[3]:
<matplotlib.text.Text at 0x7f5681ccfb10>
What is the expected correlation of parameters versus truncation mass?¶
In [4]:
fig,ax = plt.subplots(3,1,sharex=True,sharey=False,figsize=(6,7),
gridspec_kw={"hspace":0.1})
K = IdealAnalytic(log_mmin=log_mmin,logHs=hs,alpha=alpha,beta=beta,lnA=0.0)
ax[0].plot(log_mmin-hs,K.corr[:,0,1],lw=2)
ax[1].plot(log_mmin-hs,K.corr[:,1,2],lw=2)
ax[2].plot(log_mmin-hs,K.corr[:,0,2],lw=2)
#ax[0].legend(loc=0,ncol=2)
ax[0].set_ylabel(r"$\log \mathcal{H}_\star-\alpha$",fontsize=15)
ax[1].set_ylabel(r"$\alpha-\beta$",fontsize=15)
ax[2].set_ylabel(r"$\log \mathcal{H}_\star-\beta$",fontsize=15)
ax[2].set_xlabel(r"$\log_{10}\left(m_{min}/\mathcal{H}_\star\right)$",fontsize=15)
Out[4]:
<matplotlib.text.Text at 0x7f5681b27090>
Determine relationship of MRP parameters to physical ones¶
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from hmf import MassFunction
from scipy.integrate import simps
from scipy.interpolate import InterpolatedUnivariateSpline as spline
from mrpy.fitting.fit_curve import get_fit_curve
from mrpy import dndm
import os
from astropy.cosmology import Planck13
In [5]:
fig_folder=""
Our goal in this example is to find functions that relate the MRP
parameters, \((\mathcal{H}_s,\alpha,\beta)\) to the physical
parameters \((z,\Omega_m,\sigma_8)\). The results of this analysis
are already stored in the physical_dependence module of mrpy.
This example shows how those models were derived explicitly, and also
gives some ideas on how to fit theoretical curves, and the issues that
can pop up.
Figures from this example appear in MRP (2018) as figures 5,6 and 7
Our first task is to set up the fiducial cosmology and choose a fitting function. We choose the Behroozi+13 SO fitting function, with a fiducial cosmology of Planck+13 (Planck+WP 68% limits):
In [18]:
h = MassFunction(cosmo_model=Planck13, hmf_model = "Behroozi",
dlnk=0.02, lnk_max=11, lnk_min=-6.5, dlog10m=0.01)
Checking Transfer model and Resolution Parameters¶
Before we begin our actual analysis, we need to make sure of a few things concerning the models and resolution. Firstly, we want to make sure that using EH isn’t a big deal. What we do is plot the resulting mass function ratio of an EH vs CAMB model, for several redshifts (dependence on other parameters will be very small), where the mass is scaled by the log mass mode:
In [2]:
def get_mass_mode(m,dndm,scalar=True):
spl = spline(np.log10(m[dndm>0]),(m**2*dndm)[dndm>0],k=4) #k=4 necessary for getting roots next
d = spl.derivative()
if scalar:
return d.roots()[0]
else:
return d.roots()
The mass mode will not change significantly over \(\Omega_m\) or \(\sigma_8\), but will change a lot for redshift. Let’s have a look at this dependence, to make sure everything is okay:
In [19]:
fig,ax = plt.subplots(1,2,figsize=(10,5))
h.update(Mmin=3,Mmax=15)
Z = np.linspace(0,8,25)
mode = np.zeros(len(Z))
for i,z in enumerate(Z):
h.update(z=z)
mode[i] = get_mass_mode(h.m,h.dndm)
ax[0].plot(h.m,h.m**2 * h.dndm,color="k")
ax[0].set_xscale('log')
ax[1].plot(Z,mode,lw=2)
ax[1].plot(Z,13 - 0.9*Z)
ax[1].set_xlabel("Redshift",fontsize=16)
ax[0].set_xlabel("Mass",fontsize=16)
ax[1].set_ylabel("Log Mass Mode",fontsize=16)
ax[0].set_ylabel(r"$m^2 dn/dm$",fontsize=16)
Out[19]:
<matplotlib.text.Text at 0x7f2a5e564e90>
This all looks good (which means that the default resolution parameters are all good for our purposes, but we’ll check this more thoroughly soon). Furthermore, we should be able to define a more efficient function which doesn’t need to have a huge range of Mmin/Mmax:
In [3]:
def get_mass_mode_h(h,Mmin=-3,Mmax=2):
estimate = 13 - 0.9*h.z
finished = False
i = 0
while not finished and i<4:
i += 1
h.update(Mmin=estimate+1.5*Mmin*i**0.5,Mmax=estimate+1.5*Mmax*i**0.5)
mode = get_mass_mode(h.M,h.dndm,scalar=False)
if len(mode)==1 and mode[0]+Mmin > h.Mmin and mode[0]+Mmax < h.Mmax:
finished = True
else:
continue
return mode[0]
Just for testing, let’s use this function to again produce a similar plot to above.
In [21]:
def get_mass_mode_h(h):
return np.log10(h.M[np.argmax(h.m*h.dndlnm)])
Z = np.linspace(0,8,25)
mode = np.zeros(len(Z))
for i,z in enumerate(Z):
h.update(z=z)
mode[i] = get_mass_mode(h.m,h.dndm)
plt.plot(Z,mode,lw=2)
plt.xlabel("Redshift",fontsize=16)
plt.ylabel("Log Mass Mode",fontsize=16)
if fig_folder:
plt.savefig(join(fig_folder,"log_mass_mode.pdf"))
In [4]:
def recast_dndm(m,dndm,mode,Mmin=-3,Mmax=2,dm=h.dlog10m):
mvec = np.arange(mode+Mmin,mode+Mmax,dm)
sp = spline(np.log10(m[dndm>0]),np.log10(dndm[dndm>0]))
return 10**sp(mvec)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-4-b9afd81f536a> in <module>()
----> 1 def recast_dndm(m,dndm,mode,Mmin=-3,Mmax=2,dm=h.dlog10m):
2 mvec = np.arange(mode+Mmin,mode+Mmax,dm)
3 sp = spline(np.log10(m[dndm>0]),np.log10(dndm[dndm>0]))
4 return 10**sp(mvec)
NameError: name 'h' is not defined
Mmin and Mmax Determination¶
We want to choose good mass limits for each parameter set. We do this by setting a constant Mmin and Mmax w.r.t. the mass mode. By setting the mode inside these limits, we ensure a good range to define \(\beta\) in each case. Let’s have a look at the mass functions inside these limits:
In [ ]:
h.update(Mmin=-1,Mmax=18)
fig,ax = plt.subplots(2,1,sharex=True,figsize=(8,6),gridspec_kw={"hspace":0})
for z in np.linspace(0,8,9):
h.update(z=z)
mode = get_mass_mode(h.m,h.dndm)
sp = spline(np.log10(h.m[h.dndm>0]),np.log10(h.dndm[h.dndm>0]))
mvec = np.arange(mode-4,mode+4,0.02)
hvec = np.arange(-4,4,0.02)
# The following is necessary since sometimes the vectors have 1 more element
N = min(len(mvec),len(hvec))
mvec = mvec[:N]
hvec = hvec[:N]
ax[0].plot(hvec,10**sp(mvec),lw=2)
ax[1].plot(hvec,10**(2*mvec) * 10**sp(mvec),lw=2)
ax[0].set_yscale('log')
ax[0].set_ylim((1e-40,1e14))
ax[1].set_xlabel(r"$\log_{10} (m/\mathcal{H}_T)$",fontsize=20)
ax[0].set_ylabel(r"$dn/dm$",fontsize=20)
ax[1].set_ylabel(r"$m^2 \frac{dn}{dm}$",fontsize=20)
h.update(z=0)
In [ ]:
Mmin = -3
Mmax = 2
Data Generation¶
The main idea is to generate a sample of parameters, \((z,\Omega_m,\sigma_8)\), based on a realistic distribution (i.e. Planck13), then for each parameter, generate the HMF between the relevant limits. Then each HMF is to be fit with the MRP. The resulting parameters will need to be fit in Eureqa, after which we’ll port the results back here ;)
So, first up, the general parameters of the samples:
In [ ]:
# General configuration for parameter samples
N_1d = 200 # Number of samples to draw for 1D data
N_3d = 2000 # Number of samples to draw for 3D data
s8_mean = h.sigma_8
s8_sd = 0.012
Om0_mean = h.cosmo.Om0
Om0_sd = 0.017
Both \(\sigma_8\) and \(\Omega_m\) will be sampled from the Planck distribution, but \(z\) will come from a log-linear distribution (which more highly weights low redshifts, but not too much). First, define a function which will generate the samples:
In [ ]:
def get_parameter_sample(N, zbound,s8_mean=s8_mean,s8_sd=s8_sd,Om0_mean=Om0_mean,Om0_sd=Om0_sd):
sigma_8 = np.random.normal(s8_mean, s8_sd, size=N)
Om0 = np.random.normal(Om0_mean, Om0_sd, size=N)
if not hasattr(zbound, "__len__"):
z = np.repeat(zbound, N)
else:
z = np.exp(np.linspace(np.log(1 + zbound[0]), np.log(1 + zbound[1]), N)) - 1
return sigma_8, Om0, z
And a function which generates the EPS mass function from each sample.
In [ ]:
def get_hmfs(h, Om0, z, s8,Mmin=-3,Mmax=1):
# If Om0 is a scalar, then things can be sped up a bit.
try:
Om0[3]
go_for_it = False
except:
go_for_it = True
# Make sure they're all vectors
Om0 = np.atleast_1d(Om0)
z = np.atleast_1d(z)
s8 = np.atleast_1d(s8)
N = max(len(Om0),len(z),len(s8))
if len(Om0)==1: Om0 = np.repeat(Om0,N)
if len(s8)==1: s8 = np.repeat(s8,N)
if len(z)==1: z = np.repeat(z,N)
Nm = len(arange(Mmin,Mmax,h.dlog10m))
dndm = np.zeros((Nm,N))
mode = np.zeros(N)
h.update(Mmin=2,Mmax=16)
for i, (s, zz, m) in enumerate(zip(s8, z, Om0)):
if i%100==0:
print float(i) * 100 / float(len(z)), "% done"
h.update(cosmo_params={"Om0":m}, z=zz, sigma_8=s)
if go_for_it:
mode[i] = get_mass_mode(h.M,h.dndm)
else:
mode[i] = get_mass_mode_h(h)
#get vals at new mass range
dndm[:,i] = recast_dndm(h.m,h.dndm,mode[i],Mmin=Mmin,Mmax=Mmax,dm=h.dlog10m)[:Nm]
new_m = arange(mode[i]+Mmin,mode[i]+Mmax,h.dlog10m)
# for some reason, sometimes new_m is 1 bigger than needed
if len(new_m)>Nm:
new_m = new_m[:Nm]
return dndm,mode
Now, try to read the samples in from a data file, otherwise produce them again.
In [ ]:
if os.path.exists("phys_dep/raw_data.npz"):
all_data = np.load("phys_dep/raw_data.npz")
# 3d
s8_3d = all_data["s8_3d"]
Om0_3d = all_data["Om0_3d"]
z_3d = all_data["z_3d"]
dndm_3d = all_data["dndm_3d"]
Ht_3d = all_data["Ht_3d"]
# 1d
s8_1d = all_data["s8_1d"]
Om0_1d = all_data["Om0_1d"]
z_1d = all_data["z_1d"]
dndm_z_1d = all_data['dndm_z_1d']
dndm_s8_1d = all_data['dndm_s8_1d']
dndm_Om0_1d = all_data['dndm_Om0_1d']
Ht_z_1d = all_data['Ht_z_1d']
Ht_s8_1d = all_data['Ht_s8_1d']
Ht_Om0_1d = all_data['Ht_Om0_1d']
del all_data
else:
# Create Samples
s8_1d, Om0_1d, z_1d = get_parameter_sample(N_1d,(0,8))
s8_3d, Om0_3d, z_3d = get_parameter_sample(N_3d,(0,8))
# Calculate HMF's
dndm_s8_1d, Ht_s8_1d = get_hmfs(h,Om0,0,s8_1d,Mmin=Mmin,Mmax=Mmax)
dndm_Om0_1d, Ht_Om0_1d = get_hmfs(h,Om0_1d,0,s8,Mmin=Mmin,Mmax=Mmax)
dndm_z_1d, Ht_z_1d = get_hmfs(h,Om0,z_1d,s8,Mmin=Mmin,Mmax=Mmax)
print "done 1d z"
dndm_3d, Ht_3d = get_hmfs(h,Om0_3d,z_3d,s8_3d,Mmin=Mmin,Mmax=Mmax)
print "done 3d 8"
For good measure, check the distributions:
In [33]:
fig,ax = plt.subplots(1,3,figsize=(12,3))
ax[0].hist(s8_1d,normed=True,bins=20)
ax[1].hist(Om0_1d,normed=True,bins=20)
ax[0].hist(s8_3d,normed=True,bins=20)
ax[1].hist(Om0_3d,normed=True,bins=20)
ax[2].hist(z_1d,normed=True,bins=20)
ax[2].hist(z_3d,normed=True,bins=20)
ax[0].set_xlabel(r"$\sigma_8$")
ax[1].set_xlabel(r"$\Omega_m$")
ax[2].set_xlabel(r"$z$")
Out[33]:
<matplotlib.text.Text at 0x7f2a5c108090>
To me, it looks like the best choice will be something like (-3,2). This will give an upper limit of about 15.5 at z=0, which is pretty close to the largest things we expect to find.At high redshift, this will reach down to \(M_{\rm min} \approx 3\), which is incredibly small, but hey, who trusts these scales anyway?
Now save the data for use next time…
In [53]:
np.savez("phys_dep/raw_data.npz",s8_3d=s8_3d,Om0_3d=Om0_3d,z_3d=z_3d,dndm_3d=dndm_3d,Ht_3d=Ht_3d,
z_1d=z_1d,dndm_z_1d=dndm_z_1d,Ht_z_1d=Ht_z_1d,
s8_1d = s8_1d,dndm_s8_1d=dndm_s8_1d,Ht_s8_1d=Ht_s8_1d,
Om0_1d = Om0_1d,dndm_Om0_1d=dndm_Om0_1d,Ht_Om0_1d=Ht_Om0_1d)
Now it would be good to have a look at all the fits to make sure they look reasonable.
In [34]:
def plot_sample(sample,m,ax=None,ratio=False,color="k",alpha=0.3,label=None,**kwargs):
if ax is None:
ax = plt.subplot(111)
for i in range(sample.shape[1]):
if i>0:
label=None
if ratio:
ax.plot(m,sample[:,i]/sample[:,0],color=color,alpha=alpha,label=label**kwargs)
else:
ax.plot(m,sample[:,i],color=color,alpha=alpha,label=label,**kwargs)
First look at all the 1d distros
In [35]:
mvec = np.arange(Mmin,Mmax,h.dlog10m)
fig,ax=plt.subplots(3,1,sharex=True,gridspec_kw={"hspace":0.0},subplot_kw={"yscale":'log'},figsize=(8,6))
for i,sample in enumerate([dndm_s8_1d,dndm_Om0_1d,dndm_z_1d]):
plot_sample(sample,mvec,ax[i],ratio=False)
ax[2].set_xlabel(r"$\log_{10}(m/\mathcal{H}_T)$",fontsize=20)
ax[1].set_ylabel(r"$dn/dm$",fontsize=24)
ax[0].text(1,1e-13,r"$\sigma_8, z=0$",fontsize=16)
ax[1].text(1,1e-13,r"$\Omega_m, z=0$",fontsize=16)
ax[2].text(1,1e-2,r"$z$",fontsize=16)
Out[35]:
<matplotlib.text.Text at 0x7f2a3ebd5610>
These look nice and contiguous as we would expect. Now the 3d distros:
In [36]:
mvec = np.arange(Mmin,Mmax,h.dlog10m)
plot_sample(dndm_3d,mvec,ratio=False)
plt.xlabel(r"$\log_{10}(m/\mathcal{H}_T)$",fontsize=20)
plt.ylabel(r"$dn/dm$",fontsize=24)
plt.text(1,1e-13,r"$z=0$",fontsize=16)
plt.text(1,1e0,r"$z=(0,8)$",fontsize=16)
plt.yscale('log')
First we define a generic function that will fit the dndm vectors.
In [ ]:
def fit_all_mrp(sample,mode,Om0=None,Mmin=-3,Mmax=1,s=0):
par = np.zeros((sample.shape[1],4))
fit = np.zeros_like(sample)
if Om0 is None:
Om0 = np.zeros(sample.shape[1])
guesstimate = [14.5,-1.9,0.8,-40]
for i,(dndm,m) in enumerate(zip(sample.T,mode)):
M = 10**np.arange(m+Mmin,m+Mmax,h.dlog10m)[:len(dndm)]
res,curve = fit_curve(M,dndm,hs0=guesstimate[0],alpha0=guesstimate[1],
beta0=guesstimate[2],lnA0=guesstimate[3])
par[i] = res.x
fit[:,i] = curve.dndm()
guesstimate = par[i]
return par,fit
And now a couple of functions which plot the results:
In [ ]:
def plot_mrp_fit(fit_dict,ax=None,Mmin=Mmin,Mmax=Mmax,dm=h.dlog10m):
if ax is None:
ax = plt.subplot(111)
mvec = np.arange(Mmin,Mmax,dm)[:dndm_3d.shape[0]]
for i in range(dndm_3d.shape[1]):
ax.plot(mvec,fit_dict["fit_3d"][:,i]/dndm_3d[:,i],color="k",alpha=0.2)
def plot_mrp_param_dep(fit_dict):
fig,ax=plt.subplots(5,3,sharex="col",gridspec_kw={"hspace":0},figsize=(12,8))
for i in range(5):
if i<4:
ax[i,0].scatter(s8_1d,fit_dict['par_s8_1d'][:,i])
else:
ax[i,0].scatter(s8_1d,Ht_s8_1d)
ax[i,0].set_ylabel([r"$\mathcal{H}_s$",r"$\alpha$",r"$\beta$",r"$\log A$",r"$\mathcal{H}_T$"][i],fontsize=16)
for i in range(5):
if i<4:
ax[i,1].scatter(Om0_1d,fit_dict['par_Om0_1d'][:,i])
else:
ax[i,1].scatter(Om0_1d,Ht_Om0_1d)
for i in range(5):
if i<4:
ax[i,2].plot(z_1d,fit_dict['par_z_1d'][:,i])
else:
ax[i,2].plot(z_1d,Ht_z_1d)
for i in range(3):
ax[3,i].set_xlabel([r"$\sigma_8$",r"$\Omega_m$",r"$z$"][i],fontsize=16)
Here we generate the data (again, only if not found in a corresponding file already):
In [39]:
if os.path.exists('phys_dep/fits.npz'):
fits = dict(np.load("phys_dep/fits.npz"))
else:
fits = {}
fits["par_s8_1d"], fits["fit_s8_1d"] = fit_all_mrp(dndm_s8_1d , Ht_s8_1d, Mmin=Mmin,Mmax=Mmax)
fits["par_Om0_1d"],fits["fit_Om0_1d"] = fit_all_mrp(dndm_Om0_1d, Ht_Om0_1d,Mmin=Mmin,Mmax=Mmax)
fits["par_z_1d"], fits["fit_z_1d"] = fit_all_mrp(dndm_z_1d , Ht_z_1d, Mmin=Mmin,Mmax=Mmax)
fits["par_3d"], fits["fit_3d"] = fit_all_mrp(dndm_3d , Ht_3d, Mmin=Mmin,Mmax=Mmax)
np.savez("phys_dep/fits.npz",**fits)
np.savetxt("phys_dep/s8_fits.dat",np.vstack((s8_1d,fits['par_s8_1d'].T)).T)
np.savetxt("phys_dep/Om0_fits.dat",np.vstack((Om0_1d,fits['par_Om0_1d'].T)).T)
np.savetxt("phys_dep/z_fits.dat",np.vstack((z_1d,fits['par_z_1d'].T)).T)
np.savetxt("phys_dep/3d_fits.dat",np.vstack((s8_3d,Om0_3d,z_3d,fits['par_3d'].T)).T)
Now plot the resulting functions.
In [40]:
plot_mrp_fit(fits)
plot_mrp_param_dep(fits)
Return from Eureqa¶
The data we generated was input to the program “Eureqa” to determine
relationships between the physical and MRP parameters. The results were
then implemented in physical_dependence.
First, we define a helper function for putting labels above entire rows/columns of axes:
In [ ]:
def suplabel(axis,label,label_prop=None,
labelpad=5,
ha='center',va='center'):
''' Add super ylabel or xlabel to the figure
Similar to matplotlib.suptitle
axis - string: "x" or "y"
label - string
label_prop - keyword dictionary for Text
labelpad - padding from the axis (default: 5)
ha - horizontal alignment (default: "center")
va - vertical alignment (default: "center")
'''
fig = plt.gcf()
xmin = []
ymin = []
for ax in fig.axes:
xmin.append(ax.get_position().xmin)
ymin.append(ax.get_position().ymin)
xmin,ymin = min(xmin),min(ymin)
dpi = fig.dpi
if axis.lower() == "y":
rotation=90.
x = xmin-float(labelpad)/dpi
y = 0.5
elif axis.lower() == 'x':
rotation = 0.
x = 0.5
y = ymin - float(labelpad)/dpi
else:
raise Exception("Unexpected axis: x or y")
if label_prop is None:
label_prop = dict()
plt.text(x,y,label,rotation=rotation,
transform=fig.transFigure,
ha=ha,va=va,
**label_prop)
We basically run through all our samples, generating fits for them based on the parameterisations we created, then we plot each as a ratio to the true HMF. For some more clarity, we do this in several redshift bins.
In [43]:
# Import the necessary function
from mrpy.extra.physical_dependence import mrp_b13, _alpha_b13, _beta_b13, _logHs_b13, _lnA_b13
# Define each mass function based on our results
post_fit = np.zeros_like(dndm_3d)
for i,(z,m,s,mode) in enumerate(zip(z_3d,Om0_3d,s8_3d,Ht_3d)):
mvec = 10**np.arange(mode+Mmin,mode+Mmax,h.dlog10m)[:len(dndm_3d)]
post_fit[:,i] = mrp_b13(mvec,z,m,s)
# Create a figure plotting the average ratio of the functions
# to the true HMF, with shading giving the uncertainty region.
ratio = post_fit/dndm_3d
hvec = np.arange(Mmin,Mmax,0.01)[:len(mvec)]
cols = ["b","g","r","c"]
fig,ax = plt.subplots(2,2,sharey=True,sharex="col",figsize=(7,5),gridspec_kw={"hspace":0,"wspace":0},
subplot_kw={"ylim":(0.93,1.09),"xlim":(-3.1,2.2)})
for j,zrange in enumerate([(0,1),(1,3),(3,6),(6,8)]):
mask = np.logical_and(z_3d>=zrange[0], z_3d < zrange[-1])
mean = np.median(ratio[:,mask],axis=1)
qlow, qhi = np.percentile(ratio[:,mask],(16,84),axis=1)
ax[divmod(j,2)].plot(hvec,mean,lw=2,color=cols[j],label="z=(%s,%s)"%zrange)
ax[divmod(j,2)].fill_between(hvec,qlow,qhi,color=cols[j],alpha=0.3)
ax[divmod(j,2)].text(-2.5,1.04,"z=(%s,%s)"%zrange,fontsize=15)
ax[divmod(j,2)].tick_params(axis='both', which='major', labelsize=13)
if j==2:
ax[divmod(j,2)].xaxis.set_ticks([-3,-2,-1,0,1])
suplabel('x',r"$\log_{10}(m/\mathcal{H}_T)$",label_prop={"fontsize":15},labelpad=7)
suplabel("y","MRP/Behroozi",label_prop={"fontsize":15},labelpad=6)
# Save for the paper!
if fig_folder:
plt.savefig(fig_folder+"param_3d.pdf")
Also, we can re-look at the parameter dependence with included parameterisations:
In [45]:
fig,ax=plt.subplots(5,3,sharex="col",gridspec_kw={"hspace":0},figsize=(12,8))
Om0 = h.cosmo.Om0
s8 = h.sigma_8
## Dependence on s8
Hs = _logHs_b13(0,Om0,s8_1d)
alpha = _alpha_b13(0,Om0,s8_1d)
beta = _beta_b13(0,Om0,s8_1d)
lnA = _lnA_b13(0,Om0,s8_1d)
param = [Hs,alpha,beta,lnA]
for i in range(5):
if i<4:
ax[i,0].scatter(s8_1d,fits['par_s8_1d'][:,i])
ax[i,0].plot(s8_1d,param[i],color="r")
else:
ax[i,0].scatter(s8_1d,Ht_s8_1d)
ax[i,0].plot(s8_1d,Hs+np.log10((alpha+2)/beta)/beta,color="r")
ax[i,0].set_ylabel([r"$\log_{10}\mathcal{H}_s$",r"$\alpha$",r"$\beta$",r"$\ln A$",r"$\log_{10}\mathcal{H}_T$"][i],fontsize=19)
## Dependence on Om0
Hs = _logHs_b13(0,Om0_1d,s8)
alpha = _alpha_b13(0,Om0_1d,s8)
beta = _beta_b13(0,Om0_1d,s8)
lnA = _lnA_b13(0,Om0_1d,s8)
param = [Hs,alpha,beta,lnA]
for i in range(5):
if i<4:
ax[i,1].scatter(Om0_1d,fits['par_Om0_1d'][:,i])
ax[i,1].plot(Om0_1d,param[i],color="r")
else:
ax[i,1].scatter(Om0_1d,Ht_Om0_1d)
ax[i,1].plot(Om0_1d,Hs+np.log10((alpha+2)/beta)/beta,color="r")
## Dependence on z
Hs = _logHs_b13(z_1d,Om0,s8)
alpha = _alpha_b13(z_1d,Om0,s8)
beta = _beta_b13(z_1d,Om0,s8)
lnA = _lnA_b13(z_1d,Om0,s8)
param = [Hs,alpha,beta,lnA]
for i in range(5):
if i<4:
ax[i,2].scatter(z_1d,fits['par_z_1d'][:,i])
ax[i,2].plot(z_1d,param[i],color="r")
else:
ax[i,2].scatter(z_1d,Ht_z_1d)
ax[i,2].plot(z_1d,Hs+np.log10((alpha+2)/beta)/beta,color="r")
for i in range(3):
ax[4,i].set_xlabel([r"$\sigma_8$",r"$\Omega_m$",r"$z$"][i],fontsize=19)
for i in range(5):
for j in range(3):
yticks = ax[i,j].yaxis.get_major_ticks()
yticks[0].label1.set_visible(False)
yticks[-1].label1.set_visible(False)
# Save for the paper!
if fig_folder:
plt.savefig(fig_folder+"param_dep.pdf")
Use different parameterisations of the MRP and check how well they perform¶
In this example, we use the framework of the reparameterise module
to check the performance of some of the different parameterizations of
the MRP. The idea is to set up large samples of halos with the same MRP
parameters at various truncation masses, and check the covariance of the
parameters at the solution in each case, for different
parameterisations. Then compare the covariance in each case.
Figure from this example is found in MRP appendix B
Create Data¶
First we’re going to set some parameters which we’ll use throughout. These typify what we expect the HMF to be around a redshift of 0 (roughly!).
In [2]:
# Imports
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import os
from mrpy.extra import reparameterise as repar
from mrpy import MRP
from mrpy.fitting.fit_sample import SimFit
In [7]:
fig_folder=""
In [3]:
# MRP parameters
hs = 14.5
alpha = -1.85
beta = 0.72
mmins = [10.0,10.5,11.0,11.5,12.0,12.5,13.0]
V = 1
nbar = 1e6
# Number of samples
Nreal = 20
# Parameters defining outcomes
maxj = 1
In [6]:
force_recalc = False
if os.path.exists("parameterisations/data.npz") and not force_recalc:
f = np.load("parameterisations/data.npz")
gg2_mean = f['gg2_mean']
gg3_mean = f['gg3_mean']
ht_mean = f['ht_mean']
gg2_sd = f['gg2_sd']
gg3_sd = f['gg3_sd']
ht_sd = f['ht_sd']
else:
# Instantiate our final result arrays
gg2 = np.zeros((Nreal,len(mmins),3,2))
gg3 = np.zeros((Nreal,len(mmins),3,2))
ht = np.zeros((Nreal,len(mmins),3,2))
np.random.seed(1010)
for im,mmin in enumerate(mmins):
mm = MRP(None, logHs=hs, alpha=alpha,beta=beta, norm=0, log_mmin=mmin)
lnA = np.log(nbar/mm.nbar)
for i in range(Nreal):
# The following checks if there are any oddities with the covariance.
# Since the covariance is derived numerically from the analytic hessian,
# if the hessian is extremely correlated, the matrix inversion can fail.
# In such a case, we re-draw the samples.
do = True
j=0
while do and j<maxj:
logm = np.log10(mm.stats.rvs(np.random.poisson(nbar)))
#Fit the data
fit = SimFit(10**logm, mmin=10**mmin).run_downhill(hs,alpha,beta,lnA)[0].x
gg2c = repar.GG2Sample(logm=logm,log_mmin=mmin,logHs=fit[0],alpha=fit[1],beta=fit[2],lnA=fit[3])
gg2[i,im,:,0] = np.sqrt(np.diag(gg2c.cov_ratio)[:3]) * gg2c.theta_T()/gg2c.p_T()
do = np.any(np.logical_or(np.isnan(gg2[i,im,:,0]),np.isinf(gg2[i,im,:,0])))
j+=1
if j>maxj:
print "ITERATION %s ON MMIN=%s HAD NANS: "%(i,mmin)
# GG2 results
gg2[i,im,:,1] = gg2c.corr_ratio[0,1],gg2c.corr_ratio[0,2],gg2c.corr_ratio[1,2]
# GG3 results
gg3c = repar.GG3Sample(logm=logm,log_mmin=mmin,logHs=fit[0],alpha=fit[1],beta=fit[2],lnA=fit[3])
gg3[i,im,:,0] = np.sqrt(np.diag(gg3c.cov_ratio)[:3])* gg3c.theta_T()/gg3c.p_T()
gg3[i,im,:,1] = gg3c.corr_ratio[0,1],gg3c.corr_ratio[0,2],gg3c.corr_ratio[1,2]
# HT results
htc = repar.HTSample(logm=logm,log_mmin=mmin,logHs=fit[0],alpha=fit[1],beta=fit[2],lnA=fit[3])
ht[i,im,:,0] = np.sqrt(np.diag(htc.cov_ratio)[:3])* htc.theta_T()/htc.p_T()
ht[i,im,:,1] = htc.corr_ratio[0,1],htc.corr_ratio[0,2],htc.corr_ratio[1,2]
# Take the mean over the first axis
gg2_mean = np.nanmean(gg2,axis=0)
gg3_mean = np.nanmean(gg3,axis=0)
ht_mean = np.nanmean(ht,axis=0)
gg2_sd = np.nanstd(gg2,axis=0)
gg3_sd = np.nanstd(gg3,axis=0)
ht_sd = np.nanstd(ht,axis=0)
# Save the data
np.savez("parameterisations/data.npz",gg2_mean=gg2_mean,gg3_mean=gg3_mean,
ht_mean=ht_mean,gg2_sd=gg2_sd,gg3_sd=gg3_sd,ht_sd=ht_sd)
Make a Plot¶
We want to make a plot which shows the relative variance/covariance in parameters of the different forms over the different truncation masses. In the top panel we’ll plot the difference in variance of the parameters. The smaller the better. In the bottom panel, we’ll plot the correlation of the parameter combinations. In this case, the closer to 0 the better. In each case, what is plotted is the ratio of the value to that in the case of the MRP. We notate the transformation parameters as \(\mu, \nu, \delta\), and compare directly to \(\log \mathcal{H}_\star, \alpha, \beta\).
In [8]:
fig,ax=plt.subplots(2,1,sharex=True,figsize=(8,6))
ax[0].plot([np.nan],[np.nan],label="GG2",color="k",ls=":",lw=2)
ax[0].plot([np.nan],[np.nan],label="GG3",color="k",ls="--",lw=2)
ax[0].plot([np.nan],[np.nan],label="HT",color="k",lw=2)
for ip,par in enumerate([r"$\mu$",r"$\nu$",r"$\delta$"]):
col1=['r','g','b'][ip]
eb = ax[0].errorbar(mmins,gg2_mean[:,ip,0],gg2_sd[:,ip,0],ls=':',color=col1,lw=2)
eb[-1][0].set_linestyle(":")
eb=ax[0].errorbar(mmins,gg3_mean[:,ip,0],gg3_sd[:,ip,0],ls="--",color=col1,lw=2)
eb[-1][0].set_linestyle("--")
ax[0].errorbar(mmins,ht_mean[:,ip,0],ht_sd[:,ip,0],color=col1,label=par,lw=2)
for ip, cmb in enumerate([r"$\mu-\nu$",r"$\mu-\delta$",r"$\nu-\delta$"]):
col2=['k','m','c'][ip]
eb=ax[1].errorbar(mmins,gg2_mean[:,ip,1],gg2_sd[:,ip,1],ls=':',color=col2,lw=2)
eb[-1][0].set_linestyle(":")
ax[1].errorbar(mmins,gg3_mean[:,ip,1],gg3_sd[:,ip,1],ls="--",color=col2,lw=2)
eb[-1][0].set_linestyle("--")
ax[1].errorbar(mmins,ht_mean[:,ip,1],ht_sd[:,ip,1],label=cmb,color=col2,lw=2)
ax[0].set_ylim((-0.5,1.5))
ax[0].legend(loc=0,ncol=2)
ax[1].legend(loc=0,ncol=3)
plt.subplots_adjust(hspace=0.1)
ax[1].set_ylim((-2.5,1.5))
ax[-1].set_xlabel(r"$m_{\rm min}$",fontsize=18)
for i in range(2):
ax[i].tick_params(axis='both', which='major', labelsize=12)
ax[0].set_ylabel("Variance Ratio",fontsize=16)
ax[1].set_ylabel("Correlation Ratio",fontsize=16)
#Save for the paper!
if fig_folder:
plt.savefig(join(fig_folder,"reparameterisations_vs_mmin.pdf"))
Use the MRP to explore the relation between halo mass and stellar mass¶
In this example, we’ll use the MRP function for something other than pure halo mass functions. We’ll use it to find the subhalo abundance-matched (SHAM) stellar-mass halo-mass relation (i.e. how much stellar mass does a halo of mass m typically contain?).
This can be done purely numerically, but we’ll show (as per section 5 of the MRP paper) that using the analytic nature of the MRP means we can do slightly better than pure numerics.
Plots from this example appear in MRP as figures 14, 15
Motivation and numerical solution¶
The SHAM approximation for calculating the relationship between stellar and halo mass is simply to equate the cumulative number densities at every mass:
The galaxy stellar mass function (GSMF) is commonly parameterised as a double-Schechter function, which has the integral
which when equated to the MRP gives us the following equation to solve:
where we use the common definitions from the MRP paper: \(z_h = (\alpha_h + 1)/\beta\) and \(x_h = (m_h/\mathcal{H}_\star)^\beta\).
NOTE: the *subhalo* mass function is to be used here, since it is the subhalos which contain galaxies. The MRP is quite able to fit the sHMF as well.
So, we’d hope we could just go and solve this equation for \(m_\star\) as a function of \(x_h\). However, unfortunately, the solution is analytically impossible (though maybe you could try?).
The first idea then would be to solve it purely numerically, which is what we’re going to do in this section.
First, some general imports and parameters:
In [4]:
# Some imports
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from mrpy.extra.physical_dependence import mrp_params_b13
from mrpy import MRP
from scipy.interpolate import InterpolatedUnivariateSpline as spline
from scipy.optimize import newton, curve_fit
from mpmath import gammainc, log, exp
In [5]:
# Some parameters
z = 0
H0 = 67.04
mhalo = np.linspace(9,15.5,200) # vector of halo masses to plot.
Now define the MRP explicitly, based on best parameters for z=0
In [6]:
# Get relevant MRP parameters at z=0 for Planck cosmology
(logHs,alpha,beta,lnA) = mrp_params_b13(z=0)
# Create an MRP instance with those params
mrp = MRP(mhalo,logHs,alpha,beta,norm=lnA)
We’ll also want to define the GSMF explicitly. We use the data from Baldry+12 to define it:
In [7]:
# Parameters (Baldry+12 defines densities in H0=100 units and masses in H0=70.0 units)
logMs=10.66 + np.log10(0.7) + np.log10(70.0/H0)
alpha1=-0.35
log_phi1=-2.402+3*np.log10(100./H0)
alpha2=-1.47
log_phi2=-3.1+3*np.log10(100./H0)
# One way to define a double-schechter function is to add two MRP's with beta=1:
def ngtm_dblschech(m,logMs,alpha1,log_phi1,alpha2,log_phi2):
schech1 = MRP(m,logMs,alpha1,1.0,norm=np.log(10)*(log_phi1 - logMs))
schech2 = MRP(m,logMs,alpha2,1.0,norm=np.log(10)*(log_phi2 - logMs))
return schech1.ngtm() + schech2.ngtm()
# The GMSF
ngtm_baldry = ngtm_dblschech(mhalo,logMs,alpha1,log_phi1,alpha2,log_phi2)
Now that we have the left-hand and right hand sides, we can solve.
The basic idea is to find which mass in the GSMF has the same value for its \(n(>m_\star)\) as a particular mass in the MRP. You can think of this as plotting the two functions together, and then drawing horizontal arrow from one to the other (see below!).
The obvious way to do this is by using splines, and that is the best way to think about it. However, numerically, this is pretty non-robust, since then the GSMF has to be calculated over a range which covers the whole HMF, before we know what that range is.
We rather do this calculation by root-finding. So we write a function
which generally solves the double schechter function for \(x\) given
an output \(q\). Then we relate this to our actual problem. For the
root-finding, we use the gammainc function straight from mpmath,
since we need to retain as much precision as possible, to enlarge the
mass range within which this works. Also, to get more reliability, we
implement an analytic derivative.
In [8]:
def inv_inc_gamma_dbl(a1,a2,n1,n2,q,x0=0):
"""
Returns x for q = n1*Gamma(a1,x)+n2*Gamma(a2,x).
This is general, so supports a<0, but is slow, uses root finding to generate
the solution.
"""
def obj_func(x):
x=np.exp(x)
p1 = n1*gammainc(a1,x)
p2 = n2*gammainc(a2,x)
rat = float(p1/p2)
if rat>5:
func = float(log(p1) + np.log1p(1/rat) - log(q))
elif rat<0.2:
func = float(log(p2) + np.log1p(rat) - log(q))
else:
func = float(log(p1+p2) - log(q))
return func#,fprime#,fprime2
def fprime(x):
x=np.exp(x)
p1 = n1*gammainc(a1,x)
p2 = n2*gammainc(a2,x)
bottom = (p1+p2)
top = n1*x**a1 + n2*x**a2
y = exp(-x)
fprime = -float(y*top/(bottom))
return fprime
res = newton(obj_func,x0=x0,fprime=fprime)#,fprime2=True)
return np.exp(res)
def mh_to_ms_schech(log_mh,ngtm,logMs,alpha1,log_phi1,alpha2,log_phi2):
"""
Generate the SMHM relation ms(mh) for a black-box mass function and single or
double Schechter SMF.
Parameters
----------
log_mh : array_like
log10 halo masses
ngtm : array_like
Cumulative mass function.
Ms, alpha1, alpha2, log_phi1, log_phi2 : float
Double-schechter parameters. log_phi is in log10.
Returns
-------
ms : array_like
The log10 stellar masses corresponding to mh.
"""
# Make sure mh is an array
log_mh = np.atleast_1d(log_mh)
ms = np.zeros_like(log_mh)
for i,ng in enumerate(ngtm):
ms[i] = inv_inc_gamma_dbl(alpha1+1,alpha2+1,10**log_phi1,10**log_phi2,ng,x0=ms[min(i-1,0)])
return np.log10(ms)+logMs
Finally we write a convenience function that calculates the transform from mhalo to mstar in both directions, and also from either to the fraction mstar/mhalo.
In [9]:
def get_all_relations(log_mhalo, mrp, **gsmf_kwargs):
ms = mh_to_ms_schech(log_mhalo, mrp, **gsmf_kwargs)
ms_to_mh = lambda sm : spline(ms,log_mhalo)(sm)
mh_to_ms = spline(log_mhalo,ms)
ms_to_frac = spline(ms,10**(ms-log_mhalo))
mh_to_frac = spline(log_mhalo,10**(ms-log_mhalo))
return ms_to_mh, mh_to_ms,ms_to_frac,mh_to_frac
In [10]:
mstar_to_mhalo, mhalo_to_mstar, mstar_to_frac, mhalo_to_frac = get_all_relations(mhalo, mrp.ngtm(),
alpha1=alpha1,logMs=logMs,log_phi1=log_phi1,
alpha2=alpha2,log_phi2=log_phi2)
And we can make a plot. We’ll make a two-panel plot with the upper panel showing the schematic of what we’re doing, and the bottom showing the numerical result.
In [11]:
# Set up the basic figure with two subplots (one will be blank for now)
fig,ax= plt.subplots(2,1,gridspec_kw={"height_ratios":(1.6,1)})
# Plot both n(>m)
ax[0].plot(mhalo,np.log10(mrp.ngtm()),label="MRP",lw=2)
ax[0].plot(mhalo,np.log10(ngtm_baldry),label="Baldry+12 GSMF",lw=2)
# Axis stylings
ax[0].legend(loc=0)
ax[0].set_ylim((-8,1))
ax[0].set_ylabel(r"$\log_{10} n(>m)$",fontsize=15)
# Here we want to draw horizontal arrows connecting the MRP to the GSMF.
# These will signify the masses of each that have the same value of n(>m), i.e., the SMHM relation.
n1 = int(len(mhalo)*2.0/5.0)
dx1 = np.log10(mhalo_to_frac(mhalo[n1])) + 0.355
ax[0].arrow(mhalo[n1],np.log10(mrp.ngtm()[n1]),dx1,0,head_width=0.15,head_length=0.1,lw=2)
n2 = int(len(mhalo)*4.0/5.0)
dx2 = np.log10(mhalo_to_frac(mhalo[n2])) + 0.155
ax[0].arrow(mhalo[n2],np.log10(mrp.ngtm()[n2]),dx2,0,head_width=0.15,head_length=0.1,lw=2)
# Now plot the actual SMHM relation in the bottom panel
ax[1].plot(mhalo,mhalo_to_frac(mhalo),color="r",lw=2)
# Bottom axis stylings
ax[0].xaxis.set_ticks([])
ax[1].set_ylabel(r"$M_\star/M_h$",fontsize=15)
ax[1].set_xlabel(r"$m$, $\log_{10}h^{-1}M_\odot$",fontsize=15)
ax[1].set_ylim((0,0.042))
plt.subplots_adjust(hspace=0.05)
# Save for the paper!
fig.savefig("../../../mrpArticle/figures/numerical_solution.pdf")
/home/steven/miniconda3/envs/mrpy/lib/python2.7/site-packages/ipykernel/__main__.py:6: RuntimeWarning: divide by zero encountered in log10
Empirical Models from Literature¶
There are several models in the literature which attempt to parameterise the ratio \(m_star/m_h\) (i.e. to fit the red curve above). In this section, we’ll define three of those, and show how they stack up to our numerical solution.
In this case all we really need is the conversion from \(m_h\) to the fraction, \(f = m_star/m_h\). We’ll take the models of Moster+09, Mutch+13 and Behroozi+10:
In [12]:
## Moster+09 (default arguments are those given in the paper)
def mh_to_frac_moster(mh,f0=0.0282,m1=11.884,beta=1.057,gamma=0.556):
return 2*f0/(10**(-beta*(mh-m1)) +10**(gamma*(mh-m1)))
## Mutch+13
def mh_to_frac_mutch(mh,eps=0.17*0.9,mpeak=11.6,sigma=0.56):
return eps * np.exp(-(mh-mpeak)**2/sigma**2)
## Behroozi+10 -- this is a bit more complicated since Behroozi defines m_h(m_s) rather than the other way around.
def ms_to_mh_behroozi(ms,m1=10.72,m1_a=0.55,m0=12.35,m0_a=0.28,beta=0.44,
beta_a=0.18,delta=0.57,delta_a=0.17,gamma=1.56,gamma_a=2.51,
z=0):
# Get evolution of parameters
a=1./(1.+z)
m1 += m1_a*(a-1)
m0 += m0_a*(a-1)
beta += beta_a*(a-1)
delta += delta_a*(a-1)
gamma += gamma_a*(a-1)
logx = ms - m1
x = 10**logx
return m0 + beta*logx + x**delta/(1+x**-gamma) - 0.5
def mh_to_ms_func_behroozi(**kwargs):
ms = np.linspace(6,13,1000)
mh = ms_to_mh_behroozi(ms,**kwargs)
return spline(mh,ms)
def mh_to_frac_behroozi(mh,**kwargs):
msfunc = mh_to_ms_func_behroozi(**kwargs)
return 10**(msfunc(mh)-mh)
With these defined, we can fit their parameters to our numerical solution, since in this case we care more about how the model itself performs, rather than the actual values of the parameters. To do this, we just use a simple curve_fit method. We perform the fit in log-log space, since that performs the best over a large mass range.
In [13]:
# Log of the numerical solution
logf = np.log(mhalo_to_frac(mhalo))
# ---------- MOSTER --------------------------------------------------
# Define an objective function to minimize for Moster+09
def moster_obj_func(mh, *args) :
return np.log(mh_to_frac_moster(mh,*args))
# Minimize the function
moster_res = curve_fit(moster_obj_func,mhalo,logf,p0=(0.0282,11.884,1.057,0.556))[0]
# ---------- MUTCH --------------------------------------------------
# Define an objective function to minimize for Moster+09
def mutch_obj_func(mh, *args) :
return np.log(mh_to_frac_mutch(mh,*args))
# Minimize the function
mutch_res = curve_fit(mutch_obj_func,mhalo,logf,p0=(0.15,11.6,0.56))[0]
# ---------- BEHROOZI --------------------------------------------------
# Define an objective function to minimize for Moster+09
def behroozi_obj_func(mh,m0,m1,beta,delta,gamma) :
return np.log(mh_to_frac_behroozi(mh,m0=m0,m1=m1,beta=beta,delta=delta,gamma=gamma))
# Minimize the function
behroozi_res = curve_fit(behroozi_obj_func,mhalo,logf,p0=(12.35,10.72,0.44,0.57,1.56))[0]
/home/steven/miniconda3/envs/mrpy/lib/python2.7/site-packages/ipykernel/__main__.py:15: RuntimeWarning: invalid value encountered in log
Now plot the solutions against the numerical result:
In [14]:
def plotbase():
fig,ax = plt.subplots(2,1,sharex=True,gridspec_kw={"height_ratios":(2,1)},figsize=(9,6))
ax[0].plot(mhalo,mhalo_to_frac(mhalo),label="Numerical",color="k")
# Plot Moster results
ax[0].plot(mhalo,mh_to_frac_moster(mhalo,*moster_res),label="Moster+09",color="C0",lw=2)
ax[1].plot(mhalo,mh_to_frac_moster(mhalo,*moster_res)/mhalo_to_frac(mhalo),color="C0",lw=2)
# Plot Mutch results
ax[0].plot(mhalo,mh_to_frac_mutch(mhalo,*mutch_res),label="Mutch+13",color="C1",lw=2)
# Plot Behroozi Results
ax[0].plot(mhalo,np.exp(behroozi_obj_func(mhalo,*behroozi_res)),label="Behroozi+10",color="C2",lw=2)
ax[1].plot(mhalo,np.exp(behroozi_obj_func(mhalo,*behroozi_res))/mhalo_to_frac(mhalo),color="C2",lw=2)
# Plot Stylings
ax[0].legend(loc=0)
ax[0].set_yscale('log')
ax[1].set_ylim((0.4,1.2))
ax[0].set_xlim((9,16))
ax[1].grid(True)
ax[0].set_ylabel(r"$M_\star/M_h$",fontsize=18)
ax[1].set_ylabel("Fit / Numeric",fontsize=14)
ax[1].set_xlabel(r"$m$, $\log_{10}h^{-1}M_\odot$",fontsize=18)
plt.subplots_adjust(hspace=0.05)
return fig,ax
fig,ax = plotbase()
We now investigate the behaviour of the parameterisations based on what we can learn from analytic approximations to the solution. Though we cannot solve the full equation analytically for all \(m_h\), we do well to identify its behaviour in the limits.
As \(x \rightarrow 0\), we have the identity \(\Gamma(z,x) \rightarrow -x^z/z\), for \(z<0\). The speed of convergence of this limit depends heavily on the value of \(z\) (the more negative, the faster the convergence).
In this case, at small mass, the \(\alpha_2\) term dominates the \(\alpha_1\) term (which generally has positive shape parameter), so we can simply write
so that we find that the ratio is a power law:
While M13 cannot replicate this behaviour in the low-mass limit, both M09 and B10 can. Specifically, M09 is equivalent in the low-mass limit if \(\beta = p\) and \(2 f_0/m_1^p = K\), while B10 requires \(\beta = 1/(p+1)\) and \(m_1 (\sqrt{10}/m_0)^{1/\beta} = K\).
Let’s write a function to get these constants of the low-mass limit:
In [15]:
def set_lowmass_const(alpha_num,alpha_den, Ms_num,Ms_den,norm_num,norm_den,
beta_num,beta_den):
"""
Sets p, K for analytic approximation of low-mass relationship.
Parameters
----------
alpha_num, alpha_den : float
slopes of the numerator and denominator mass functions respectively
(eg. for SMHM, this would alpha_smf, and alpha_hmf).
Ms_num, Ms_den : float
Characteristic log10 masses of numerator and denominator mf's.
(eg. for SMHM, this would be Ms, Hs)
norm_num, norm_den : float
Normalisation of numerator and denominator mf's.
(eg. for SMHM, this would be phi, A*10^Hs)
beta_num, beta_den : float
Cut-off parameter for num and den mf's.
(eg. for SMHM, this would be 1, beta)
"""
z_num = (alpha_num + 1)/beta_num
z_den = (alpha_den + 1)/beta_den
p = (alpha_den-alpha_num)/z_num
scale = 10**(Ms_num-Ms_den)
z_rat = z_num/z_den
N = norm_den/norm_num
K = scale * (N * z_rat)**(1.0/z_num)*10**(-Ms_den*p)
return p,K
In [16]:
# Actually set the p,K in our case
p,K = set_lowmass_const(alpha2,alpha, logMs,logHs,10**log_phi2,np.exp(lnA)*10**logHs,1.0,beta)
Now we can re-write our objective functions for Moster+09 and Behroozi+10 so that their constants are set by the low-mass behaviour, and rederive solutions.
In [17]:
def moster_fixed_obj_func(mh,m1,gamma):
f0 = 10**(m1*p) * K/2
beta = p
return np.log(mh_to_frac_moster(mh,f0,m1,beta,gamma))
moster_fixed_res = curve_fit(moster_fixed_obj_func,mhalo,logf,p0=(11.884,0.6))[0]
def behroozi_fixed_obj_func(mh,m0,delta,gamma):
beta = 1./(p+1)
m1 = np.log10(K) + (m0-0.5)/beta
return np.log(mh_to_frac_behroozi(mh,m1=m1,m0=m0,beta=beta,delta=delta,gamma=gamma))
behroozi_fixed_res = curve_fit(behroozi_fixed_obj_func,mhalo,logf,p0=(12.35,0.57,1.56))[0]
Plot the solutions along with the plot from before. For some reason in the notebook the figures get closed after every cell so we have to recall the earlier plot.
In [23]:
fig,ax = plotbase()
def fixed_additions():
ax[0].plot(mhalo,np.exp(moster_fixed_obj_func(mhalo,*moster_fixed_res)),color="C0",ls='--',lw=2)
ax[1].plot(mhalo,np.exp(moster_fixed_obj_func(mhalo,*moster_fixed_res))/mhalo_to_frac(mhalo),color="C0",ls="--",lw=2)
ax[0].plot(mhalo,np.exp(behroozi_fixed_obj_func(mhalo,*behroozi_fixed_res)),color="C2",ls="--",lw=2)
ax[1].plot(mhalo,np.exp(behroozi_fixed_obj_func(mhalo,*behroozi_fixed_res))/mhalo_to_frac(mhalo),color="C2",ls="--",lw=2)
# add label to fig
ax[0].plot([0],[0],color='k',ls='--',label="Fixed low-mass")
return fig, ax
fig,ax = fixed_additions()
We can see clearly that the asymptotic behaviour at low mass for the dashed curves we just added tends towards unity.
Defining our own extended model¶
The most severe deficiency of both models, and this is exacerbated by correctly setting the low-mass behaviour, is the position of the turning point. A simple way to change the position of the turning point without affecting the behaviour in the limits, is to use an extension of M09:
where \(w\) controls the peak position, \(k\) is able to correct the amplitude of the high-mass power-law, and \(\delta\) adds the flexibility needed to induce an upturn left of the turning point. This extension obeys the same relations as M09 in terms of the low-mass approximation, but has the difference that the high-mass power law has the slope \(\delta-\gamma\).
We can define this function, and its objective function for optimization:
In [18]:
def mh_to_frac_triple_pl(mh,f0=0.0282,m1=11.884,beta=1.057,gamma=0.556,delta=0,k=0,w=1):
return (10**(w*delta*(mh-m1)) +2*f0)/(10**(-beta*(mh-m1)) +w*10**(gamma*(mh-m1))+k)
def triple_pl_obj_func( mh,m1,gamma,delta,k,w):
f0 = 10**(m1*p) * K/2
beta = p
return np.log(mh_to_frac_triple_pl(mhalo,f0,m1,beta,gamma,delta,k,w))
triple_pl_res = curve_fit(triple_pl_obj_func,mhalo,logf,p0=(11.884,2.6,2,1,1))
In [21]:
fig,ax = plotbase()
#fig,ax = fixed_additions()
ax[0].plot(mhalo,np.exp(triple_pl_obj_func(mhalo,*triple_pl_res[0])),label="TPL",color="C3",ls='--',lw=2)
ax[1].plot(mhalo,np.exp(triple_pl_obj_func(mhalo,*triple_pl_res[0]))/mhalo_to_frac(mhalo),color="C3",ls="--",lw=2)
ax[0].legend(loc=0,ncol=2)
ax[0].set_ylim((1e-4,0.06))
ax[0].set_xlim((9,15.5))
ax[1].set_ylim(0.8,1.2)
# Save for the paper!
fig.savefig("../../../mrpArticle/figures/compare_mine_fixed.pdf")
Construct a custom fit for an extension of the MRP: a double-Schechter function.¶
Note: this example is almost identical to that found in the R package ``tggd`` in the ``tggd_log`` documentation.
While mrpy has several in-built routines which aid in fitting the
MRP function to data, it does not natively support fitting extensions
of the MRP, such as double-Schechter functions. For such functions, one
can fairly simply create custom fits using the methods found in Scipy,
for example.
In this example, we create a double-Schechter galaxy stellar mass function (GSMF) down to a target stellar mass (xmin) of log10(SM) = 8. We use data from Baldry+2012 to define the function:
Both mixtures have
Mixture 1 has:
Mixture 2 has:
Furthermore, we use only the purely statistical routines of mrpy to
achieve our results:
In [1]:
# Imports
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gaussian_kde
from scipy.optimize import minimize
from mrpy.base import stats
First we define objects which capture the statistical quantities of each mixture:
In [2]:
mix1 = stats.TGGDlog(10.66,-1.47,1.0,8)
mix2 = stats.TGGDlog(10.66,-0.35,1.0,8)
\(\phi^\star\) is defined as the value of the pdf in log-space at the pivot scale by \(e\), i.e.
This normalisation is important in our sampling since it gives the ratio of samples produced by each mixture. We can produce the relevant normalisation using:
In [3]:
M1norm=0.79/mix1.pdf(10.66)
M2norm=3.96/mix2.pdf(10.66)
Mtot=M1norm+M2norm
Now say we would like to sample 1e5 galaxies, we can produce these like so:
In [4]:
Nsamp=1e5
np.random.seed(100)
mix1_sample = mix1.rvs(int(Nsamp*M1norm/Mtot))
mix2_sample = mix2.rvs(int(Nsamp*M2norm/Mtot))
gal_sample = np.concatenate((mix1_sample,mix2_sample))
Let’s plot the distribution we’ve just created, to check if everything has worked okay:
In [5]:
# Create x array to plot against
x = np.linspace(8,12,401)
plt.plot(x,gaussian_kde(gal_sample)(x),label="Sample", lw=2)
plt.plot(x,mix1.pdf(x) * M1norm/Mtot,label="Mixture 1", lw=2)
plt.plot(x,mix2.pdf(x) * M2norm/Mtot,label="Mixture 2", lw=2)
plt.plot(x,mix1.pdf(x) * M1norm/Mtot + mix2.pdf(x) * M2norm/Mtot, label="GSMF", lw=2)
plt.yscale('log')
plt.ylim((1e-7,1))
plt.legend(loc=0)
Out[5]:
<matplotlib.legend.Legend at 0x7fa87a9d5210>
Now we can try to fit the mixed model. The trick here is we fit for the
mixture using an additional parameter \(\lambda\), where one
component is multiplied by \(\lambda\) and the other
\(1-\lambda\). We define it so that M1norm/Mtot=lambda and
M1norm/Mtot=1-lambda. Here’s our model:
In [6]:
def mix_lnl(par,data):
mix1 = stats.TGGDlog(par[0],par[1],1.0,8)
mix2 = stats.TGGDlog(par[0],par[2],1.0,8)
return -np.sum(np.log(mix1.pdf(data)*par[3] + mix2.pdf(data)*(1-par[3])))
Now we can perform the fit, using a downhill-gradient method of our choice. The fit is probably not fantastic though. Generalised Gamma distributions (including truncated ones) display poor convergence properties using ML. Full MCMC is a better route when trying to fit GSMF type data. And the data certainly should not be binned!
In [7]:
GSMFfit = minimize(mix_lnl, x0=[10,-2,0,0.5], args=(gal_sample,), bounds=[(9,12),(-2.5,-0.5),(-1.0,0.5),(0.0,1.0)])
print "Maximum likelihood parameters: ", GSMFfit.x
Maximum likelihood parameters: [ 10.65427592 -1.46844631 -0.34935003 0.83445411]
This accords very well with our input parameters!
License¶
The MIT License (MIT)
Copyright (c) 2016 Steven Murray
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Changelog¶
Development Version¶
v1.1.0 [8th Jan 2018]¶
This version is the version used for all plots in Murray, Robotham, Power (2018), and is released along with that paper. There are many changes in the code from previous versions, the result of a couple of years of sporadic work.
v1.0.0¶
Features¶
- New
PerObjFitclass supersedesget_fit_perobjfunction, providing more coherent fitting capabilities. - Added heaps of “real-world” examples (used in MRP paper):
- https://github/steven-murray/mrpy/docs/examples/fit_curve_against_analytic.ipynb
- https://github/steven-murray/mrpy/docs/examples/fit_simulation_suite.ipynb
- https://github/steven-murray/mrpy/docs/examples/heirarchical_model_stan.ipynb
- https://github/steven-murray/mrpy/docs/examples/explore_analytic_model.ipynb
- https://github/steven-murray/mrpy/docs/examples/mmin_dependence.ipynb
- https://github/steven-murray/mrpy/docs/examples/physical_dependence.ipynb
- https://github/steven-murray/mrpy/docs/examples/parameterization_performance.ipynb
- https://github/steven-murray/mrpy/docs/examples/SMHM.ipynb
- Added
modelargument tofit_perobj_stanto facilitate pickling of multiple fits. - Added ability to send keyword arguments to priors in
PerObjFitclass - Added a
normal_priorfunction for simple normal priors.
Enhancements¶
- Changed default weighting from 1 to 0 in
get_fit_curve. - Added tests for the
PerObjLikeWeightsclass. - Added tests for
nbarandrhobarfor generalmin ``MRP` subclasses. - Changed imports so that they wouldn’t show up in docs
- Many improvements to documentation (including this file!)
Bugfixes¶
- Fixed issue setting
log_mmininIdealAnalytic - Fixed issue in which
nbarandrhobarare wrong ifmminis notm.min()inMRPsubclasses.
API Summary¶
mrpy.base.special |
Definitions of all special functions used throughout mrpy. |
mrpy.base.stats |
A module defining the TGGD distribution (as well as in log and ln space) in standard R style. |
mrpy.base.core |
Basic MRP functionality, such as functions to generate the MRP with different normalisations. |
mrpy.extra.physical_dependence |
Module containing functions for the dependence of MRP parameters on physical parameters, defined with respect to Behroozi+13. |
mrpy.extra.likelihoods |
Provides classes which extend the basic mrpy.core.MRP class. |
mrpy.extra.analytic_model |
A module defining the likelihoods and related quantities involved when the data is purely ideal and analytic. |
mrpy.extra.reparameterise |
Provides classes which implement variations of the MRP, in which the parameters have been transformed. |
mrpy.fitting.fit_curve |
Routines that implement simple least-squares fits directly to dn/dm without errors. |
mrpy.fitting.fit_sample |
Routines that implement fits directly to samples of masses (or other values drawn from a TGGD). |