Fitting Mass Functions with emcee

In this demo, we provide a very simple example of how to fit mass functions with emcee. This example does not do anything with the complexity required for fitting real-world data. It does not account for measurement error or Eddington bias. However, it gives a taste of what is achievable.

[87]:

import emcee
import hmf
import numpy as np

import matplotlib.pyplot as plt
from scipy.stats import poisson
from multiprocessing import Pool



from IPython.display import display, Markdown

%matplotlib inline

[2]:

emcee.__version__

[2]:

'3.0.2'

[3]:

hmf.__version__

[3]:

'3.1.3.dev4+ge4212c6.d20200721'

Create Some Mock Data

First, let’s create a default model mass function:

[43]:

model = hmf.MassFunction(z=1.5, transfer_model='EH')

Now, let’s create some mock data, based on the true mass function, but with some randomness assuming the mass function was measured (perfectly) in some volume:

[44]:

volume = 100 ** 3 # Mpc^3 / h^3

[45]:

counts_per_bin = np.random.poisson(model.dndlog10m * volume * model.dlog10m)

[46]:

plt.plot(model.m, model.dndlog10m * volume * model.dlog10m)
plt.plot(model.m, counts_per_bin)
plt.xscale('log')
plt.yscale('log')
../_images/examples_fitting_11_0.png

Define a likelihood

Now let’s define a likelihood function, based on some input model for dndlog10m. The model is simply a Poisson model. We take the model into the space in which it describes an actual number of observations, and get the log pdf there:

[47]:

def log_likelihood(model, data, volume, dlog10m):
    return np.sum(poisson.logpmf(data, mu=model*volume*dlog10m))

Define an emcee-compatible likelihood function

Now we define a likelihood function for emcee. There is a bit more flexibility here, as this function needs to calculate priors on all input parameters, and handle exceptions as well. This means this function is rather specific to the problem at hand. We will define a very simple function, but one that is fairly general.

[48]:

fiducial_model = model.clone()

First, we define a small utility function that will take a dictionary in which keys may be dot-paths, and converts it to a nested dictionary:

[49]:

def flat_to_nested_dict(dct: dict) -> dict:
    """Convert a dct of key: value pairs into a nested dict.

    Keys that have dots in them indicate nested structure.
    """
    def key_to_dct(key, val, dct):
        if '.' in key:
            key, parts = key.split('.', maxsplit=1)

            if key not in dct:
                dct[key] = {}

            key_to_dct(parts, val, dct[key])
        else:
            dct[key] = val

    out = {}
    for k, v in dct.items():
        key_to_dct(k, v, out)

    return out

So, this will do the following:

[50]:

flat_to_nested_dict(
    {
        'nested.key': 1,
        'nested.key2': 2,
        'non_nested': 3

    }
)

[50]:

{'nested': {'key': 1, 'key2': 2}, 'non_nested': 3}

This will enable us to pass a list of parameter names that we want updated, which could be parameters of nested models. This means our posterior function is fairly general, and can accept any model parameters to be updated:

[98]:

def log_prob(param_values, param_names, data, model, volume, derived=None):
    # Update the base model with all the parameters that are being constrained.
    params = flat_to_nested_dict(dict(zip(param_names, param_values)))

    model.update(**params)

    if derived is not None:
        derived = [getattr(model, d) for d in derived]
    else:
        derived = []

    return log_likelihood(model.dndlog10m, data, volume, model.dlog10m), derived

We can test that the log_prob function works:

[52]:

log_prob([0.9], ['sigma_8'], counts_per_bin, model, volume, derived=['mass_nonlinear'])

[52]:

(-1390.6879386397318, [array(5.83204315e+10)])

Notice the derived parameters: we can pass any quantity of the MassFunction, and it will be stored on every iteration. Nice!

Run emcee

Let’s first run a simple model in which we just want to fit \(\sigma_8\) and the spectral index \(n_s\). We use the popular emcee package, and pass in our log_prob model:

[60]:

sampler = emcee.EnsembleSampler(
    nwalkers = 100,
    ndim = 2,
    log_prob_fn = log_prob,
    kwargs = {
        'param_names': ['sigma_8', 'n'],
        'data': counts_per_bin,
        'model': model,
        'volume': volume,
        'derived': ['mass_nonlinear']
    },
    pool = Pool(4),
)

/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:377: UserWarning: Nonlinear mass outside mass range
  warnings.warn("Nonlinear mass outside mass range")
/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:377: UserWarning: Nonlinear mass outside mass range
  warnings.warn("Nonlinear mass outside mass range")
/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:377: UserWarning: Nonlinear mass outside mass range
  warnings.warn("Nonlinear mass outside mass range")

On the advice of the emcee documentation, we set up some initial positions of the walkers around the solution.

[61]:

initialpos = np.array([fiducial_model.sigma_8, fiducial_model.n]) + 1e-4 * np.random.normal(size=(sampler.nwalkers, sampler.ndim))

[62]:

sampler.run_mcmc(initialpos, nsteps=1000)

[62]:

State([[0.812517   0.96740204]
 [0.82984608 0.96694073]
 [0.81653361 0.9656326 ]
 [0.82261185 0.9665801 ]
 [0.8102518  0.96749993]
 [0.82020062 0.96884229]
 [0.81968681 0.96613499]
 [0.80751869 0.96396036]
 [0.81132914 0.96522728]
 [0.81345488 0.96482091]
 [0.81815306 0.96822774]
 [0.82297313 0.96462598]
 [0.82597001 0.96222415]
 [0.8204828  0.96716804]
 [0.82306268 0.96922072]
 [0.82241968 0.96456152]
 [0.82463299 0.9655332 ]
 [0.82045326 0.96579144]
 [0.81639143 0.96658896]
 [0.82310428 0.9682027 ]
 [0.82192434 0.96464978]
 [0.83147352 0.96654571]
 [0.8157654  0.96777175]
 [0.83319911 0.96809932]
 [0.81681679 0.9635151 ]
 [0.82759967 0.96613953]
 [0.82531001 0.96399171]
 [0.82192877 0.96582809]
 [0.81542795 0.9673587 ]
 [0.81230205 0.96544857]
 [0.81980445 0.96677717]
 [0.82473735 0.9661013 ]
 [0.83556881 0.96525546]
 [0.8145986  0.96269262]
 [0.81476969 0.96651966]
 [0.82230772 0.96429339]
 [0.82630066 0.96775328]
 [0.8170044  0.9620626 ]
 [0.81458535 0.96246315]
 [0.81471739 0.96628634]
 [0.82734829 0.96353397]
 [0.82123505 0.96576662]
 [0.81452712 0.96796586]
 [0.81843941 0.96565461]
 [0.81678365 0.96191043]
 [0.82880826 0.96322063]
 [0.82690041 0.96280437]
 [0.81234069 0.96687445]
 [0.81833558 0.96742033]
 [0.81188034 0.96350916]
 [0.82064682 0.96507912]
 [0.83007665 0.9632581 ]
 [0.82343703 0.96941583]
 [0.82139026 0.96612144]
 [0.82500533 0.96379157]
 [0.81868078 0.96902344]
 [0.82051868 0.96682397]
 [0.82243481 0.96584191]
 [0.82013074 0.96382395]
 [0.81619164 0.96515733]
 [0.81984934 0.96751159]
 [0.81569339 0.96636135]
 [0.83262378 0.96858608]
 [0.82537676 0.96507569]
 [0.82169244 0.96399202]
 [0.82259506 0.96954585]
 [0.81162461 0.96412451]
 [0.81190128 0.9665054 ]
 [0.8175165  0.96546521]
 [0.81006349 0.96652593]
 [0.81088687 0.96368845]
 [0.8238553  0.96701842]
 [0.82488918 0.96619036]
 [0.81048038 0.96460308]
 [0.82614518 0.9660271 ]
 [0.8235058  0.96689262]
 [0.81577499 0.96506154]
 [0.82445992 0.96511738]
 [0.82241107 0.96609058]
 [0.81213763 0.9660414 ]
 [0.81028155 0.9654621 ]
 [0.81940399 0.97021029]
 [0.82693119 0.96653072]
 [0.81862029 0.9668193 ]
 [0.82067371 0.9647356 ]
 [0.82839192 0.963874  ]
 [0.809108   0.96696672]
 [0.82132415 0.96543191]
 [0.82574774 0.96393896]
 [0.82616052 0.96638275]
 [0.81915099 0.96582024]
 [0.82237358 0.96515534]
 [0.83688351 0.96418875]
 [0.81852464 0.96585457]
 [0.81285049 0.96711565]
 [0.8199895  0.9647657 ]
 [0.81963206 0.96338418]
 [0.81767408 0.96475152]
 [0.81698525 0.96946209]
 [0.81725398 0.96785498]], log_prob=[-1321.6468555  -1322.52336774 -1320.82115653 -1321.01426989
 -1322.13347053 -1322.35533506 -1320.74680145 -1323.53951807
 -1321.68404161 -1321.39797316 -1321.72784215 -1321.09999745
 -1323.26345382 -1321.08899881 -1323.03532484 -1321.07741124
 -1321.11800667 -1320.74276535 -1320.90964786 -1322.00366832
 -1321.01013801 -1322.8495541  -1321.47842399 -1324.33998589
 -1321.74275509 -1321.70434785 -1321.66743377 -1320.82510487
 -1321.26810679 -1321.41964293 -1320.9010349  -1321.16769159
 -1324.1709755  -1322.8413994  -1321.0379521  -1321.18429957
 -1322.15985558 -1323.269938   -1323.11384998 -1321.00627943
 -1322.30422001 -1320.77800489 -1321.71384606 -1320.73004281
 -1323.50045722 -1322.85079113 -1322.81290782 -1321.47379661
 -1321.17279566 -1322.58079765 -1320.82331432 -1323.14018705
 -1323.3221903  -1320.81442667 -1321.74093293 -1322.49509933
 -1320.94415875 -1320.86695212 -1321.36823306 -1320.92795258
 -1321.24785381 -1320.91722817 -1324.58610898 -1321.28006987
 -1321.30690166 -1323.38225282 -1322.15660864 -1321.47195975
 -1320.77869767 -1321.8789143  -1322.68042721 -1321.3250135
 -1321.20598589 -1322.15580027 -1321.39141283 -1321.22540074
 -1320.99385637 -1321.13799991 -1320.88815269 -1321.39023989
 -1321.87715249 -1324.0565578  -1321.65890442 -1320.9013828
 -1320.92133778 -1322.31341947 -1322.21457395 -1320.7965959
 -1321.7641621  -1321.46525572 -1320.71979122 -1320.9062473
 -1325.00247426 -1320.7251741  -1321.47033812 -1320.89659695
 -1321.695903   -1320.94860038 -1323.00891678 -1321.46047235], blobs=[2.43237797e+10 2.91648801e+10 2.49626622e+10 2.69081466e+10
 2.37519275e+10 2.68154812e+10 2.59597826e+10 2.22294542e+10
 2.34865733e+10 2.39415325e+10 2.60720089e+10 2.64934764e+10
 2.67136068e+10 2.64539805e+10 2.77522888e+10 2.63193424e+10
 2.72132211e+10 2.60856815e+10 2.51646358e+10 2.74876658e+10
 2.62024084e+10 2.95559695e+10 2.52923527e+10 3.05543612e+10
 2.45110797e+10 2.82562169e+10 2.69955015e+10 2.65129271e+10
 2.50958840e+10 2.37920233e+10 2.61596404e+10 2.73971170e+10
 3.04713151e+10 2.37252439e+10 2.47078485e+10 2.62174546e+10
 2.83140956e+10 2.42029199e+10 2.36664933e+10 2.46358528e+10
 2.74656153e+10 2.62997588e+10 2.50033912e+10 2.54895362e+10
 2.41073681e+10 2.78104202e+10 2.71374291e+10 2.41480616e+10
 2.59141286e+10 2.32195362e+10 2.59550348e+10 2.81988036e+10
 2.79166915e+10 2.64368908e+10 2.68540629e+10 2.64287198e+10
 2.63738402e+10 2.66611057e+10 2.54885107e+10 2.47513868e+10
 2.63642129e+10 2.49176218e+10 3.05140807e+10 2.73066039e+10
 2.59659188e+10 2.77027067e+10 2.33000050e+10 2.39427699e+10
 2.51880449e+10 2.34699281e+10 2.30088789e+10 2.73871390e+10
 2.74657502e+10 2.31208894e+10 2.77914560e+10 2.72511091e+10
 2.46156049e+10 2.70512011e+10 2.67202090e+10 2.38923119e+10
 2.32731549e+10 2.69495059e+10 2.81642520e+10 2.58385227e+10
 2.58736187e+10 2.78665609e+10 2.33289425e+10 2.62373743e+10
 2.71081525e+10 2.78935801e+10 2.57290381e+10 2.64623221e+10
 3.05758631e+10 2.55641719e+10 2.43420792e+10 2.56908462e+10
 2.52395103e+10 2.50509288e+10 2.60639090e+10 2.57244997e+10], random_state=('MT19937', array([3010991687, 1433366889, 4221405966, 1416151253, 3571036671,
       3571567351, 3617448558, 1708944356, 2355426697, 1975937285,
       3274553095, 3023766973, 4077549661, 3789930581,  873399172,
        868365483, 3296907483,  334032045, 3808706278, 1441354804,
       2014058982, 3607909725, 3452910208, 1411274757, 2062595218,
       2876891164, 1385668807, 4052497944,  665092010, 2423041690,
        395260579, 1282554622,  870191133, 3933196700, 1396143140,
       1095897011, 2369267201,  657731204,  834372367, 3284579739,
        219773595, 1330355438, 3375518106, 2711778082, 1837927017,
       1204990710, 4042827915, 1069892769, 2242376285,  661762990,
       3898970883, 3092654721,  906073352, 1255118679,  677245811,
       3825484840, 1924460280,  246112358, 1003177144, 1400168998,
       3502415521,  917281527, 4084550239, 1919567022, 2069619878,
        741611488, 3503439977,  193452069, 4104958293, 4247830910,
       3829228084, 1864232335, 2347051935, 1437298699, 2519712600,
       3064761557, 2590550356, 1743410357, 3760085702, 1709383524,
       1751293765, 1473007862, 1852997295, 3912295010, 3989041097,
       1927332117, 1001638882, 3684246907, 1206769616, 1276227959,
       4106449237, 3040437620, 4068953156, 2591667847, 4174197124,
        809558309,  770698902,  639505468, 1136190391, 1635456729,
        779995844, 3154881102, 2301101302, 2414368414, 3313271264,
       1581666314, 1387140706, 3482159655, 1547374700, 1831709399,
       2980407292, 1438636775, 2157209452, 1586032405, 3390334548,
        993219402, 1113086213,  709758385,  748490091, 3745675161,
        138706538, 3460914348, 3169812733, 2256278506,  623177898,
       3017175548, 3562368808,  904548435, 1110064479,  389050167,
        252102857, 1563989322, 1180302627,  459055110,  627989599,
       3018748411, 1451339223,  959784469, 3097074112, 2514658713,
       2464622003, 2087079152, 1867651254, 2717430748, 3299737911,
       3249936435,  825233657,  455511082, 2823889225, 4113124419,
       3089413335, 3578929134,  423837805,  150075681, 3772034753,
       3420921347, 2761688900, 1099324841, 2837475302, 3743183947,
       4077423891,  488799992, 2412862165, 3966439414, 3439170851,
         81188597,   49497071, 3840244400, 4040237744, 2891235300,
          3442727, 2836977957, 3486096401, 2238833608, 3487809619,
       3155152280, 4113000771, 3508727274, 3958557665, 1548861139,
       3402192982, 2233272614, 2195152730,  178735003, 1186388937,
       3959849889,  749318738,  392727489, 1106483805, 2556275863,
       3033016435, 3740444587, 1109541441, 3437238003,  909491960,
       3556188170,  189073105, 1121202647, 1072427671,  927260018,
       1824592484, 2750270138, 2573358151, 1885091053,  928025271,
        631173641,  845485739,  977127828, 1431784584,  858557931,
       4146320444, 3358534732, 3110633641, 1444678819,  254417518,
       4084754866, 4140526091, 2056544203,  216004678, 1242755478,
        453583462, 3718330882,  517236880, 1700461219, 3942453166,
       2648019520, 1776667548, 3313827839, 2301162440,  886357793,
        603316446, 1467257211, 2998702311,  210510884, 1042598417,
       3429242311, 2960440244, 1283222150, 1924889872, 2993775698,
       2471074342,  821293982, 1716000998, 2691036463, 1327531324,
        593573618, 3978099819, 4277429419, 2078607857,  593918363,
       4292036669, 3313375763, 2411311131, 3897805207, 2257206579,
       4069444456, 2742366635, 2516455334, 1130445333, 3344296704,
       2013760045,   98183799,  569022296, 3665793712, 2803882733,
        781782741,  185465298, 3929878580, 4234698322,   87745905,
       3755386722, 2308050559, 2793944942, 3148749364,  581001756,
       1547662942, 3416250046, 1339357844, 2396838141, 4272094892,
       1128502434, 3900946728, 3311338200, 3307175279, 1640891606,
       4100456435, 1745950520, 3967814615, 3532941775, 3082305130,
       2359549487, 4207761885, 2606420013,  102797824, 3055493477,
        394259133, 2039453341, 1877336789, 2340998722,  780799928,
        938814388, 2771433303, 3196454979, 3352703017, 1126407398,
       1530383391, 1092098769, 2714043020,  421011266, 2536176501,
       2226686724,  245473205, 3171003372, 3443182864,  599676995,
       1096705436,  971412639,  675815764, 2638098655,  253441355,
       1271857756, 3460681547, 2043527436, 1101169785, 3639565441,
       3553622061, 2910599094, 2782404150,   71011899, 2185989828,
       2059799318, 2574732472, 3838606186,  604894924,  314022656,
       4259692150, 2249781825, 3414642279, 3199016856, 3295823582,
       4267407365,  723611958, 2578396183, 1573397202, 3239680964,
         85026493, 4159867595, 1127793160, 1034939957, 3115819481,
        665096369, 1369370115, 1141236003, 3728346089, 1299368033,
        649332664, 3377083031, 2416029668,  303403121, 2310508304,
       2683484722, 2275657629, 3088922205, 2992749925, 3326099467,
       2888032499, 2891724994, 3313575824, 4167813906, 3256413369,
       1187504249,  369615155, 2439934992, 3388295164, 3668652842,
       2017075645, 2254628167, 3028960353,  917230102, 3589371289,
       2205353917, 4094540016, 2537826809, 3501831714,  136333403,
       4171261624, 1512837059, 2653939944, 3265001615, 3973707302,
       1137124443, 2887283226,  305745721,  703413495, 1230917737,
       3447653671,  637571428, 1938827366,  421456484, 1639748825,
       1091489420, 1350083345, 3468074647, 1798931265, 3437590154,
       3338560785, 3164670085,  825293874, 1280690567, 4194750462,
       3557118509, 3091135836, 3027604326, 2880680358,  938049006,
       4280523127, 3252425678, 1104263493,   81338268, 3606248064,
       2886333570, 2067133002, 1580435540, 1345254171,  695882776,
       3316965053, 2571705372, 3418551749,  827938975, 1295436624,
         76296670, 2295760466, 3303053892, 1462904270, 4153691213,
       3708139390, 2623634494,  892247839, 2795798641, 2669384127,
       4017377120,   74316956, 4085698396, 2910077659,  412899045,
       2964087952, 1843263204, 3481911643, 4172483608, 2392633725,
       2637039959, 1038236072, 2686420500, 3884299265,  682801512,
        828362509,  782282506, 2616171383, 3418205957, 2094639530,
       4021917379,  354313191, 2180560934, 4115179089, 1561190080,
       2620981921,  365432673, 2357114351,  843484284, 1702920698,
       1111277404, 1327463776, 2397409494, 1804785312, 4055920856,
       1352761305, 1933930167,   31723046, 2772243848, 4269071588,
       3411620890, 1475275535, 2663099780,   87819342, 2178690480,
       2103465561, 1306577459, 1194113820, 3881869916, 4078664281,
       4238588340, 3115971636, 3938695065,  553473943, 3550106407,
       2479809016,   20931116,  429327547,  818525329, 1284250903,
        113601256,  819638460, 2196423883,   23207807,  975749962,
       2089670681, 1109917774,   58482065,  169920096, 3152479287,
        127530720, 2660207730,   41164155, 3305460700, 2785121729,
       1862495723,  227213191, 2831217540, 2926377220, 2135577429,
       1589696410, 2488923588, 1636149543, 3323017083, 4165709967,
       2874083036, 1889135079, 2607084371, 1789986751, 2380325660,
        670719651, 1224293234, 3691150884,   41449699, 2267935622,
       3919979405, 1315635034,   88504060, 4149622353, 1849789313,
        917831583, 1379938533, 1201309631, 1916489413,  939972893,
       2559681176, 2608651579,  128801680, 3336439102, 4205359464,
       3828736694, 2019738972,  615246111,  229646520, 1698302061,
       1842255150,  274408184, 3649941914,  730146959, 1893410131,
        407308887, 2690056582, 2652547863,  872367158, 1492485614,
       2167076284, 1350370028, 2006375714, 3950254470, 2771767598,
       1318091799, 2469195389, 1979294351, 1365918688, 1183647078,
        566655925, 1183675245,  586625897,  243677354, 2389983177,
       3241818314, 1161339473,  480728130, 2829007285,  184202686,
       2648706587, 2429610140, 2878526130, 1118671801,  927419511,
       3137099824, 1109566077,   19619536, 3011767601, 1515942761,
        978363680, 2258141828, 1147637173, 3035709180, 1793921134,
       2471664108,  709617420, 1679336914, 3530543110, 1969798978,
       1713618494, 3740529277,    7127194, 3596808460, 3909109339,
       2198411294, 4154633975, 1492379234, 2701612560, 1431699733,
       3276341236, 4076082523, 1795101047, 3788458884, 3655229412,
         20755759, 2809107836, 3574786795, 2796825240], dtype=uint32), 345, 0, 0.0))

Now we can plot the posterior as a heatmap, with the true input values shown as white lines:

[78]:

plt.hexbin(sampler.chain[:,500:, 0].flatten(), sampler.chain[:,500:, 1].flatten())
plt.axvline(fiducial_model.sigma_8, color='white')
plt.axhline(fiducial_model.n, color='white')
plt.colorbar()

[78]:

<matplotlib.colorbar.Colorbar at 0x7f90c68e6d30>
../_images/examples_fitting_34_1.png

And we’re done! The posterior contains the truth to within 1-sigma, and we have ~1% error-bars on \(\sigma_8\) and ~0.2% error-bars on \(n\).

An example of fitting a new fitting function

Perhaps you are a simulator and want to deliver a new fitting function, or at least new parameters for a particular fitting function. Let’s look at an example of this. We choose the Sheth-Tormen function:

[94]:

fiducial_model = hmf.MassFunction(z=0, hmf_model='SMT', dlog10m=0.2, Mmin=10, Mmax=16)

The functional form of \(f(\nu)\) is:

[91]:

display(Markdown(f"${fiducial_model.hmf._eq}$"))

\(A\sqrt{2a/\pi}\nu\exp(-a\nu^2/2)(1+(a\nu^2)^{-p})\)

The fitted parameters are:

[92]:

fiducial_model.hmf.params

[92]:

{'a': 0.707, 'p': 0.3, 'A': 0.3222}

Again, we use a Poisson model to create mock data, this time in a Gpc volume:

[95]:

volume = 500 ** 3 # Mpc^3 / h^3

counts_per_bin = np.random.poisson(fiducial_model.dndlog10m * volume * fiducial_model.dlog10m)

plt.plot(fiducial_model.m, fiducial_model.dndlog10m * volume * fiducial_model.dlog10m)
plt.plot(fiducial_model.m, counts_per_bin)
plt.xscale('log')
plt.yscale('log')
../_images/examples_fitting_44_0.png

We’ve already set up our probability function, which was already general enough to be able to fit the hmf parameters, so off we go!

[96]:

smt_model = fiducial_model.clone()

[99]:

sampler_hmf = emcee.EnsembleSampler(
    nwalkers = 50,
    ndim = 3,
    log_prob_fn = log_prob,
    kwargs = {
        'param_names': ['hmf_params.a', 'hmf_params.p', 'hmf_params.A'],
        'data': counts_per_bin,
        'model': smt_model,
        'volume': volume,
        'derived': ['dndlog10m']
    },
    pool = Pool(4),
)

[100]:

init_pos = (
    np.array([
        smt_model.hmf.params['a'],
        smt_model.hmf.params['p'],
        smt_model.hmf.params['A']]) +
    1e-4 * np.random.normal(size=(sampler_hmf.nwalkers, sampler_hmf.ndim))
)

[103]:

sampler_hmf.run_mcmc(init_pos, nsteps=300, progress='notebook')



[103]:

State([[0.70571872 0.30004471 0.32231539]
 [0.70126824 0.30036688 0.32263641]
 [0.70542706 0.30004243 0.32235664]
 [0.70898865 0.30002223 0.32194784]
 [0.70514584 0.30010418 0.32234736]
 [0.70760904 0.30015674 0.32201378]
 [0.70865575 0.29998153 0.32200619]
 [0.70629783 0.29995717 0.32231429]
 [0.70670178 0.30012446 0.32214238]
 [0.70399317 0.30004328 0.32251483]
 [0.71196335 0.29985674 0.32170366]
 [0.71287314 0.29977231 0.32165765]
 [0.70657605 0.30000788 0.322248  ]
 [0.70154182 0.30031702 0.32263111]
 [0.70613473 0.30017031 0.32218223]
 [0.69897046 0.30028188 0.32297026]
 [0.71286506 0.2998973  0.3215743 ]
 [0.70397023 0.30011033 0.322486  ]
 [0.70337736 0.30024625 0.32245751]
 [0.70244604 0.30028499 0.32254351]
 [0.70931581 0.30014018 0.32182071]
 [0.7067235  0.29986426 0.32232561]
 [0.70624881 0.3001988  0.3221465 ]
 [0.70516839 0.3000907  0.32235771]
 [0.70520897 0.30030826 0.32219668]
 [0.7095473  0.30007336 0.3218437 ]
 [0.70706776 0.30020669 0.32204248]
 [0.70424787 0.30013549 0.32243317]
 [0.7067557  0.30012488 0.322138  ]
 [0.70295151 0.30030533 0.32245804]
 [0.7038777  0.30004237 0.32253163]
 [0.70854828 0.2998389  0.32212853]
 [0.7093244  0.3000901  0.32186687]
 [0.70970097 0.30007077 0.32182694]
 [0.70529725 0.3001928  0.32226819]
 [0.71016046 0.2999964  0.3218182 ]
 [0.70443266 0.30011627 0.32240762]
 [0.7106056  0.29994202 0.32180073]
 [0.70851975 0.30002726 0.32198968]
 [0.70434435 0.30031496 0.32229049]
 [0.7084999  0.30003641 0.32199076]
 [0.70149862 0.30023359 0.3226938 ]
 [0.70933192 0.30008167 0.32185246]
 [0.70729214 0.30005568 0.32211522]
 [0.71175351 0.2999096  0.32169285]
 [0.70682625 0.30014886 0.32210615]
 [0.70661983 0.3003228  0.32201294]
 [0.70439053 0.30010849 0.32243311]
 [0.70949498 0.29998618 0.32190394]
 [0.70629854 0.3000359  0.32226316]], log_prob=[-291.97328516 -295.94346599 -292.31931005 -291.74488602 -291.73765271
 -291.80725566 -292.22600405 -293.84609132 -291.38487631 -294.52721767
 -293.83724984 -295.66411856 -292.9420158  -293.73660519 -291.5098835
 -296.84738261 -293.76350197 -292.71065455 -292.29132272 -292.95507237
 -292.83729568 -296.16663894 -291.8550809  -292.40072465 -293.63463677
 -291.90615767 -292.40546252 -291.99716017 -291.27751233 -294.42569631
 -294.11456616 -295.73798005 -293.14878539 -291.97085974 -291.75317897
 -292.33082701 -294.54427609 -293.15184193 -292.09940775 -293.69266687
 -291.45281226 -293.36738657 -293.44319877 -292.07758313 -293.07549881
 -292.26442973 -295.60072648 -292.00028676 -292.17693171 -293.02857205], blobs=[[4.36155563e+01 2.84344925e+01 1.85272590e+01 ... 8.38507034e-08
  9.56372526e-09 5.35470813e-10]
 [4.36192349e+01 2.84360206e+01 1.85277228e+01 ... 8.68583571e-08
  1.00235459e-08 5.70407680e-10]
 [4.36161706e+01 2.84349038e+01 1.85275360e+01 ... 8.40514526e-08
  9.59399627e-09 5.37738285e-10]
 ...
 [4.36163696e+01 2.84348571e+01 1.85273986e+01 ... 8.47453091e-08
  9.69963734e-09 5.45720476e-10]
 [4.36150607e+01 2.84342836e+01 1.85271710e+01 ... 8.13304485e-08
  9.18441092e-09 5.07189625e-10]
 [4.36170030e+01 2.84354528e+01 1.85278918e+01 ... 8.34616107e-08
  9.50481121e-09 5.31046486e-10]], random_state=('MT19937', array([2695209238, 2647629147, 2193354619, 2649505322, 3246661341,
         51232968,   79962949,  179873701, 3983381240, 3370185750,
        445723221,  346033715, 3274079981, 2602548037, 1957702839,
       3461532576, 4161222576, 1038082278,  878450754, 1855409407,
        724066930, 2311148021, 3891813334, 3263803156, 2853706301,
       2847586980, 1671173016,  735939060,  429185192, 1363064350,
       2138848637,  850664440, 1792314194,  430557638, 2733706743,
       1207343644, 3496929848, 3881131080,  392740563, 4124936328,
        456456251, 3530026234, 2439348864, 3869244314,  377897866,
        770961074, 1213315404,  195611045,  254775657,  170437308,
       1266045128, 3112887178, 4081760005, 1669808052,  439342391,
       2449766217, 2759665021, 1621364076, 2707099658, 2804233581,
         18189609,  411693864, 3953138838, 2310649322, 3447542007,
       2287034944, 2025620103, 2250466890, 2389997835,  247527463,
       2082249112, 2820535931,  388301994,  305278516, 1264971166,
        110868511,  834774623, 2143469080, 2157183775,  706972368,
       2600695753, 3990090487,   29908668, 2348379960, 3909521533,
       3272447149, 3652130046, 1342559821, 1031067751, 3621803399,
       2779121456, 3352784358, 2042722360,  122920656, 2268213812,
       3482702161, 3438593671, 2318455723, 2455999334, 2138268658,
       1425420258, 2502229046, 1064394732, 2570400070,  746152391,
       3624890135, 1519459636,  638555637, 1389789579,  747484829,
       2475391687, 1498611689, 1003994038, 3341165672, 1670614520,
       3081647314,   92426843, 1453750730, 3550983803,  898249345,
       4269689740, 2674049461, 4275288442, 3153561100,  136691932,
       3340127591,  121861724,  947633261, 1456776912, 4125252795,
       3559683688,  557446984,  452533039, 1497278549, 2891972626,
       2836336200,  304000778, 1292339171,  946954854, 1338650136,
       1180859292, 1549344673, 3145802359, 2521125079, 2137220248,
       1640630604, 3423736792,  343934670, 1392888990, 4204591814,
       4210917169, 3342206687, 3772986540, 3521113294, 1473243150,
        644386694,  469127328, 1534736211, 3447892343, 2905226407,
       2652836610, 1275310234,  597817947, 3963881130,  662742347,
       4288860228,   41194792, 4174082572,  447169300, 1265630298,
       3066861970, 3996788221,  608183650, 3024185488, 2369500720,
       1090041449, 3939304911,  925037035, 1765086068, 1337744520,
       2866256580,  270294378, 4108713331, 2296602491,  254951702,
       2384066128, 3638385403,     712310, 1882208272,  192794138,
       3773442483, 3244057607, 1948353659, 2055456869,  900881485,
       1631301931,   25050028, 3052277837, 1402012007, 2397867662,
        816280189,  818764046, 3509355338,  606661196, 2857408573,
        242422348,  924196866,   59786504, 2790724941, 2659379008,
       1205222488, 1618503389, 3136805518, 3136195442, 1868190649,
       2888368716, 3434006518, 1190434689, 2906490750, 3431509262,
       1074642057, 4038132221, 3651550126, 2376139706, 3797613524,
       1227896495, 1997599976, 2468315774,  279657761,  844741445,
        617827114, 1099964320, 3667159088,  537018706, 2631468182,
       2168024562, 2299460868,  822525943, 3458927839,  864761243,
       2232687879, 2018528050,  922555956, 3152383575, 4137296225,
        433672236, 1414486742,  302349133, 2615420069, 2665160146,
       1072914862, 1087699620,   21346057, 1121602451, 1610964641,
       1092251061,  626972149, 2514265657, 1897926980, 4261144573,
       4060741681, 1542945183, 2649030300, 3527857413, 4098343921,
       1254376557, 2684640795, 2108481555, 2165763767, 3129138705,
       2124008782, 2923461739, 1983134045,  849013318, 1688004886,
        907512555, 2907583218, 1791350905,   41927223, 1264631153,
       4278809323, 3841931267, 1033782221, 2995494388,    4755582,
       3701567992, 3739998142,  695732873, 1084890294, 3034717033,
       3626648469,  695507250, 2269468095,  194485190, 2802248601,
       2648171266, 2566903326,  691052568, 3863138291, 2739798308,
        566444817,  400987765, 3922583635, 1948668692, 1995863976,
       3583671899, 3333747608, 1704741624, 1593207028, 4259776371,
       4133861978, 2905210027, 2500477908, 3917348634, 1032203945,
       1655178619, 3688631023, 1703645483, 3449857533, 3301588642,
        422497768,  872574394, 3457381080, 3228852806, 2514371337,
       1042588975, 4020226419, 1893455269, 3591940824, 4279950885,
       4178011202, 3405262077, 1145025427, 3513649778, 3794908252,
       1284594804, 1172203724, 3860875240, 2664775555, 2823845950,
       1340050203, 3627077542, 3539951354,  344802387, 2194421369,
       3281407244, 1995346482, 1079071201, 2807374371, 1598396607,
        261568444,     326839, 2825325092, 4090002079, 4047939121,
       3826124547, 2477010580, 3380211954, 2727965387, 3234635433,
        549733799, 2649691353,  366209251, 4228503155, 1987874040,
       2474903336, 1889126683, 2151236132, 3726858731,  425076785,
       2968544695, 2075851894, 2831996298,   63581575,  345516043,
       4193784166, 1767562833, 3728434816, 1199931048, 3115746379,
       3676241015, 3173393321, 3295024441, 2364485778, 2366428860,
       1275547801, 2354495069,   83730810,  718383532, 2786935316,
       2719873924, 2907524249,  136110385, 2717158293,  796096497,
       3923001143, 2254244844, 2991923610, 4086704724, 3725676723,
        307291181, 1200373680, 1620627900, 1380452724, 1222848420,
        917498110, 3415519546, 2757001008, 2127283752,   86884430,
        299884896, 4139500276, 3137696915, 2115850591, 1147152690,
        743165082, 1545887490, 1258100288, 4288035638, 4196988703,
       4171528593, 1183856867, 1455609119, 3470687830, 2225919801,
        787272135,  279320834, 4036072850,  304889672, 1726999843,
        307605912,  616898197,  801113494, 3210307683,   49615787,
       3460255429, 1655112827, 4026304912, 1770208335,  934437779,
       2855949828,  267008067, 3298131547, 3034016266, 2432814301,
       2912640166, 4272465227,   24182798, 4294545685,   64849515,
       2818315199,  601186570,   82836877, 2407288428, 2715774129,
       3890918037, 2321031121, 2627730700, 2796361440,   65202128,
       2719193199, 3985400800, 3570056396, 1825150553, 4030777940,
       4287015212, 1606870516, 2025127433, 3072872497,  464572871,
        114080185, 3169698464, 2720260534, 2223198208, 3272695669,
       3471164299, 1501740301, 1728792602, 2734157041, 2127893324,
       1786695425, 1618439030, 1516331467, 1002298222, 2410399620,
       4285972770,  142229308, 3945572832, 3900998899, 2794080412,
       4159512065, 1414727068, 3301714505,  144334919,  265912629,
       3530354009, 2024671445, 2717774424, 4065126470, 1243127162,
        963165084, 2493127522, 1758543534, 2619984954, 3400615260,
       3322033734,  635724031,  294449927,  812640229,  849031458,
       4087285948, 2859677201, 1231886253, 3263696041, 2020295525,
       4129796440, 2166724510, 3038989539, 1990376252,   87243479,
       2881630872, 3857552180, 1277212483, 1209268234, 4188376748,
       3352406736, 2672146901, 1208060266, 3766116369,  986437198,
       1442155787,  347929042, 2711347629,  934880859, 2499788704,
       1464802702, 3502557472, 3126249404, 3368939045, 1167039125,
       2344332556, 3475424878, 2131621501, 3105809766, 2040152474,
       2859813022, 4059534384,  340083373,  103878588, 3115311293,
       2205194488, 1084077185, 4149393776, 1192712839, 2988805805,
       3565986741, 3722273086, 2153314963, 2884624353, 2580215900,
       3899518847,  829130844, 3982076591, 3147994997, 2804004889,
       2717379546, 4050846669, 1597829971,  834448887, 1859330038,
       2881001046, 2746277752, 3457348486,   46229092, 3073076409,
        353130301, 2548238153, 2202603845, 4087219187, 2379793796,
       1400039462,  574217339, 3745870590, 2568676665, 4204914693,
       4178420999, 1025049129, 2312079962, 1529366137,  807512890,
       3313511292, 3016205122, 3458450123, 3872882951, 2889320531,
        849298377, 3777060733, 1968905376,  178914076, 1285520316,
        506607041, 2752741930, 2752335761, 1852914679, 3789250202,
       2354055305,  566936168, 4246993117,  755219712, 2082541872,
       3551269501, 3035511209, 4008947182, 3838599958, 3825540922,
        906345092, 1860802238, 4218213997, 4107261563, 4199247801,
       1288820235, 1129275429, 4200462492, 3591363431], dtype=uint32), 600, 0, 0.0))

Again, we can plot the posterior triangle plot:

[109]:

fig, ax = plt.subplots(2, 2, sharex='col', sharey='row', figsize=(10, 10), gridspec_kw={"hspace":0, 'wspace': 0})

flatchain = sampler_hmf.chain[:, 100:, :].reshape((-1, 3))

ax[0, 0].hexbin(flatchain[:, 0], flatchain[:, 1])
ax[1, 0].hexbin(flatchain[:, 0], flatchain[:, 2])
ax[1, 1].hexbin(flatchain[:, 1], flatchain[:, 2])

ax[1, 0].set_xlabel("a", fontsize=15)
ax[1, 0].set_ylabel("A", fontsize=15)
ax[1, 1].set_xlabel("p", fontsize=15)
ax[0, 0].set_ylabel("p", fontsize=15)


[109]:

Text(0, 0.5, 'p')
../_images/examples_fitting_51_1.png

Since we collected \(dn/d\log_{10}m\) as a derived parameter, we can also plot it:

[116]:

dndlog10m_samples = sampler_hmf.blobs[100:, :, :].reshape((-1, 40))

We’ll make a plot of the posterior distribution of the mass function, as a fraction of the median of the distribution:

[118]:

percentiles = np.percentile(dndlog10m_samples, q=(16, 50, 84), axis=0)

[128]:

plt.fill_between(fiducial_model.m, percentiles[0]/percentiles[1], percentiles[2]/percentiles[1])
plt.xscale('log')
plt.ylabel("Fractional difference to Median")
plt.xlabel("Mass")
#plt.yscale('log')

[128]:

Text(0.5, 0, 'Mass')
../_images/examples_fitting_56_1.png