Model fitting 4: Only Synchrotron

import warnings
warnings.filterwarnings('ignore')

import matplotlib.pylab as plt
import jetset
from jetset.test_data_helper import  test_SEDs
from jetset.data_loader import ObsData,Data
from jetset.plot_sedfit import PlotSED
from jetset.test_data_helper import  test_SEDs

print(jetset.__version__)

1.2.2

test_SEDs

['/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_3C345.ecsv',
 '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk421_EBL_DEABS.ecsv',
 '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv',
 '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_DEABS.ecsv']

Loading data

see the data_format user guide for further information about loading data

print(test_SEDs[1])
data=Data.from_file(test_SEDs[1])

/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk421_EBL_DEABS.ecsv

%matplotlib inline
sed_data=ObsData(data_table=data)
sed_data.group_data(bin_width=0.2)

sed_data.add_systematics(0.1,[10.**6,10.**29])
p=sed_data.plot_sed()

================================================================================

*  binning data  *
---> N bins= 89
---> bin_widht= 0.2
================================================================================
../../../_images/Jet_example_only_synchrotron_8_1.png 
sed_data.save('Mrk_401.pkl')

Phenomenological model constraining

see the Phenomenological model constraining: application user guide for further information about loading data

Spectral indices

from jetset.sed_shaper import  SEDShape
my_shape=SEDShape(sed_data)
my_shape.eval_indices(silent=True)
p=my_shape.plot_indices()
p.setlim(y_min=1E-15,y_max=1E-6)

================================================================================

* evaluating spectral indices for data *
================================================================================
../../../_images/Jet_example_only_synchrotron_13_1.png

Sed shaper

mm,best_fit=my_shape.sync_fit(check_host_gal_template=False,
                  Ep_start=None,
                  minimizer='lsb',
                  silent=True,
                  fit_range=[10., 21.])

================================================================================

* Log-Polynomial fitting of the synchrotron component *
---> first blind fit run,  fit range: [10.0, 21.0]
---> class:  HSP
Table length=4
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-1.545300e-01-1.545300e-019.534795e-03---1.000000e+00-1.000000e+010.000000e+00False
LogCubicc-1.023245e-02-1.023245e-021.433073e-03---1.000000e+00-1.000000e+011.000000e+01False
LogCubicEp1.672267e+011.672267e+014.139942e-02--1.667039e+010.000000e+003.000000e+01False
LogCubicSp-9.491659e+00-9.491659e+002.515285e-02---1.000000e+01-3.000000e+010.000000e+00False
---> sync       nu_p=+1.672267e+01 (err=+4.139942e-02)  nuFnu_p=-9.491659e+00 (err=+2.515285e-02) curv.=-1.545300e-01 (err=+9.534795e-03)
================================================================================

my_shape.IC_fit(fit_range=[23., 29.],minimizer='minuit',silent=True)
p=my_shape.plot_shape_fit()
p.setlim(y_min=1E-15)

================================================================================

* Log-Polynomial fitting of the IC component *
---> fit range: [23.0, 29.0]
---> LogCubic fit
Table length=4
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-2.098186e-01-2.098186e-013.133032e-02---1.000000e+00-1.000000e+010.000000e+00False
LogCubicc-4.661868e-02-4.661868e-022.178352e-02---1.000000e+00-1.000000e+011.000000e+01False
LogCubicEp2.524926e+012.524926e+011.147759e-01--2.529412e+010.000000e+003.000000e+01False
LogCubicSp-1.011085e+01-1.011085e+013.498963e-02---1.000000e+01-3.000000e+010.000000e+00False
---> IC         nu_p=+2.524926e+01 (err=+1.147759e-01)  nuFnu_p=-1.011085e+01 (err=+3.498963e-02) curv.=-2.098186e-01 (err=+3.133032e-02)
================================================================================
../../../_images/Jet_example_only_synchrotron_16_3.png

Model constraining

from jetset.obs_constrain import ObsConstrain
from jetset.model_manager import  FitModel
from jetset.minimizer import fit_SED
sed_obspar=ObsConstrain(beaming=25,
                        B_range=[0.001,0.1],
                        distr_e='lppl',
                        t_var_sec=3*86400,
                        nu_cut_IR=1E12,
                        SEDShape=my_shape)


prefit_jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=False,silent=True)
prefit_jet.save_model('prefit_jet.pkl')

================================================================================

*  constrains parameters from observable *
Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm3.112712e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss5.050000e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.080000e-020.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.697542e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.373160e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm39.060842e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.188500e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.181578e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature7.726502e-01-1.500000e+011.500000e+01FalseFalse
================================================================================

pl=prefit_jet.plot_model(sed_data=sed_data)
pl.add_residual_plot(prefit_jet,sed_data)
pl.setlim(y_min=1E-15,x_min=1E7,x_max=1E29)
../../../_images/Jet_example_only_synchrotron_19_0.png

Model fitting only Synchorotron component

Note

With the new implementation of composite model (FitModel class) to set parameters you have to specify the model component, this is different from versions<1.1.2, and this holds also for the freeze method and for setting fit_range intervals, and for the methods relate to parameters setting in general. See the Composite Models and depending pars user guide for further information about the new implementation of FitModel, in particular for parameter setting

Model fitting with Minuit

from jetset.jet_model import Jet
jet_minuit=Jet.load_model('prefit_jet.pkl')
jet_minuit.set_gamma_grid_size(200)
Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.697542e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.373160e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm39.060842e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.188500e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.181578e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature7.726502e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm3.112712e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss5.050000e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.080000e-020.000000e+00--FalseFalse

we switch off the IC component

jet_minuit.spectral_components.SSC.state='off'
jet_minuit.show_model()

--------------------------------------------------------------------------------
model description:
--------------------------------------------------------------------------------
type: Jet
name: jet_leptonic

electrons distribution:
 type: lppl
 gamma energy grid size:  201
 gmin grid : 4.697542e+02
 gmax grid : 1.373160e+06
 normalization:  True
 log-values:  False
 ratio of cold protons to relativistic electrons: 1.000000e-01

radiative fields:
 seed photons grid size:  100
 IC emission grid size:  100
 source emissivity lower bound :  1.000000e-120
 spectral components:
   name:Sum, state: on
   name:Sync, state: self-abs
   name:SSC, state: off
external fields transformation method: blob

SED info:
 nu grid size jetkernel: 1000
 nu size: 500
 nu mix (Hz): 1.000000e+06
 nu max (Hz): 1.000000e+30

flux plot lower bound   :  1.000000e-30

--------------------------------------------------------------------------------
Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.697542e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.373160e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm39.060842e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.188500e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.181578e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature7.726502e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm3.112712e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss5.050000e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.080000e-020.000000e+00--FalseFalse
--------------------------------------------------------------------------------

fit_model_minuit=FitModel( jet=jet_minuit, name='Only-Synch-best-fit-minuit')

fit_model_minuit.freeze('jet_leptonic','z_cosm')
fit_model_minuit.freeze('jet_leptonic','R_H')
fit_model_minuit.freeze('jet_leptonic','R')
fit_model_minuit.freeze('jet_leptonic','gmax')
fit_model_minuit.jet_leptonic.parameters.R.fit_range=[10**15.5, 10**17.5]
fit_model_minuit.jet_leptonic.parameters.beam_obj.fit_range=[5., 50.]

from jetset.minimizer import fit_SED,ModelMinimizer

model_minimizer_minuit=ModelMinimizer('minuit')
best_fit_minuit=model_minimizer_minuit.fit(fit_model_minuit,sed_data,10.0**15,10**20.0,fitname='SSC-best-fit-minuit',repeat=3)

filtering data in fit range = [1.000000e+15,1.000000e+20]
data length 13
================================================================================

* start fit process *
-----
fit run: 0
0it [00:00, ?it/s]

- best chisq=1.36478e+00

fit run: 1
- old chisq=1.36478e+00

0it [00:00, ?it/s]

- best chisq=1.36477e+00

fit run: 2
- old chisq=1.36477e+00

0it [00:00, ?it/s]

- best chisq=1.36477e+00

-------------------------------------------------------------------------
Fit report

Model: SSC-best-fit-minuit
Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*9.157430e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.373160e+061.000000e+001.000000e+15FalseTrue
jet_leptonicNemitters_density1 / cm31.850602e+000.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*4.329106e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.159322e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature6.029154e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm3.112712e+161.000000e+031.000000e+30FalseTrue
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss2.090652e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*2.326797e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.080000e-020.000000e+00--FalseTrue
converged=True
calls=95
mesg=
FCN = 1.365 Nfcn = 95
EDM = 2.08e-05 (Goal: 0.0002)
Valid Minimum Valid Parameters No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok APPROXIMATE NOT pos. def. FORCED
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 par_0 0.92e3 0.27e3 1 1E+09
1 par_1 1.9 0.5 0
2 par_2 43e3 8e3 1 1E+09
3 par_3 2.16 0.06 -10 10
4 par_4 0.60 0.08 -15 15
5 par_5 0.0209 0.0023 0
6 par_6 23.3 1.6 5 50
par_0 par_1 par_2 par_3 par_4 par_5 par_6
par_0 7.19e+04 -87.1 (-0.604) 9.91e+04 (0.044) 4.55 (0.297) -4.59 (-0.228) -0.0167 (-0.028) -174 (-0.398)
par_1 -87.1 (-0.604) 0.289 90.4 (0.020) 0.00427 (0.139) -0.00639 (-0.159) -0.000193 (-0.160) -0.105 (-0.120)
par_2 9.91e+04 (0.044) 90.4 (0.020) 7.15e+07 224 (0.464) 427 (0.674) 0.75 (0.039) 294 (0.021)
par_3 4.55 (0.297) 0.00427 (0.139) 224 (0.464) 0.00327 -0.000393 (-0.092) 7.76e-06 (0.060) 0.00549 (0.059)
par_4 -4.59 (-0.228) -0.00639 (-0.159) 427 (0.674) -0.000393 (-0.092) 0.00562 6.13e-05 (0.364) 0.0139 (0.114)
par_5 -0.0167 (-0.028) -0.000193 (-0.160) 0.75 (0.039) 7.76e-06 (0.060) 6.13e-05 (0.364) 5.06e-06 -0.00146 (-0.397)
par_6 -174 (-0.398) -0.105 (-0.120) 294 (0.021) 0.00549 (0.059) 0.0139 (0.114) -0.00146 (-0.397) 2.66 
dof=6
chisq=1.364771, chisq/red=0.227462 null hypothesis sig=0.967966

best fit pars
Table length=12
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
jet_leptonicgmin9.157430e+029.157430e+022.681484e+02--4.697542e+021.000000e+001.000000e+09False
jet_leptonicgmax1.373160e+06------1.373160e+061.000000e+001.000000e+15True
jet_leptonicN1.850602e+001.850602e+005.363195e-01--9.060842e-010.000000e+00--False
jet_leptonicgamma0_log_parab4.329106e+044.329106e+048.457976e+03--3.188500e+041.000000e+001.000000e+09False
jet_leptonics2.159322e+002.159322e+005.717871e-02--2.181578e+00-1.000000e+011.000000e+01False
jet_leptonicr6.029154e-016.029154e-017.496457e-02--7.726502e-01-1.500000e+011.500000e+01False
jet_leptonicR3.112712e+16------3.112712e+163.162278e+153.162278e+17True
jet_leptonicR_H1.000000e+17------1.000000e+170.000000e+00--True
jet_leptonicB2.090652e-022.090652e-022.250326e-03--5.050000e-020.000000e+00--False
jet_leptonicNH_cold_to_rel_e1.000000e-01------1.000000e-010.000000e+00--True
jet_leptonicbeam_obj2.326797e+012.326797e+011.628276e+00--2.500000e+015.000000e+005.000000e+01False
jet_leptonicz_cosm3.080000e-02------3.080000e-020.000000e+00--True
-------------------------------------------------------------------------

================================================================================

%matplotlib inline
fit_model_minuit.set_nu_grid(1E6,1E30,200)
fit_model_minuit.eval()
p2=fit_model_minuit.plot_model(sed_data=sed_data)
p2.setlim(y_min=1E-11,x_min=1E15,y_max=1E-9,x_max=3E19)
../../../_images/Jet_example_only_synchrotron_28_0.png 
try:
    c=model_minimizer_minuit.minimizer.draw_contour('r','s')
except:
    pass

m=model_minimizer_minuit.minimizer.draw_profile('r')
../../../_images/Jet_example_only_synchrotron_30_0.png 
best_fit_minuit.save_report('SSC-best-fit-minuit.pkl')
model_minimizer_minuit.save_model('model_minimizer_minuit.pkl')
fit_model_minuit.save_model('fit_model_minuit.pkl')

MCMC sampling

from jetset.mcmc import McmcSampler
from jetset.minimizer import ModelMinimizer

We used a flat prior centered on the best fit value. Setting bound=5.0 and bound_rel=True means that:

  1. the prior interval will be defined as [best_fit_val - delta_m , best_fit_val + delta_p]

  2. with delta_p=delta_m=best_fit_val*bound

If we set bound_rel=False then delta_p = delta_m = best_fit_err*bound

It is possible to define asymmetric boundaries e.g. bound=[2.0,5.0] meaning that

  1. for bound_rel=True

    delta_p = best_fit_val*bound[1]

    delta_m =b est_fit_val*bound[0]

  2. for bound_rel=False

    delta_p = best_fit_err*bound[1]

    delta_m = best_fit_err*bound[0]

In the next release a more flexible prior interface will be added, including different type of priors

Given the large parameter space, we select a sub sample of parameters using the use_labels_dict. If we do not pass the ‘use_labels_dict’ the full set of free parameters will be used

model_minimizer_minuit = ModelMinimizer.load_model('model_minimizer_minuit.pkl')

mcmc=McmcSampler(model_minimizer_minuit)

labels=['N','B','beam_obj','s','gamma0_log_parab']
model_name='jet_leptonic'
use_labels_dict={model_name:labels}

mcmc.run_sampler(nwalkers=128,burnin=10,steps=50,bound=5.0,bound_rel=True,threads=None,walker_start_bound=0.005,use_labels_dict=use_labels_dict)

mcmc run starting

0%|          | 0/50 [00:00<?, ?it/s]

mcmc run done, with 1 threads took 232.73 seconds

print(mcmc.acceptance_fraction)

0.52171875

p=mcmc.plot_model(sed_data=sed_data,fit_range=[1E15, 1E20],size=100)
p.setlim(y_min=1E-14,x_min=1E6,x_max=2E20)
../../../_images/Jet_example_only_synchrotron_37_0.png 
p=mcmc.plot_model(sed_data=sed_data,fit_range=[1E15, 1E20],size=100,quantiles=[0.05,0.95])
p.setlim(y_min=1E-14,x_min=1E6,x_max=2E20)
../../../_images/Jet_example_only_synchrotron_38_0.png

f=mcmc.plot_chain(‘s’,log_plot=False)

f=mcmc.corner_plot()
../../../_images/Jet_example_only_synchrotron_40_0.png 
mcmc.get_par('N')

(array([1.82210412, 1.84695012, 1.80917268, ..., 1.91776914, 1.76293129,
        1.91006206]),
 0)

f=mcmc.plot_par('beam_obj')
../../../_images/Jet_example_only_synchrotron_42_0.png 
f=mcmc.plot_par('gamma0_log_parab',log_plot=True)
../../../_images/Jet_example_only_synchrotron_43_0.png

Save and reuse MCMC

mcmc.save('mcmc_sampler.pkl')

from jetset.mcmc import McmcSampler
from jetset.data_loader import ObsData
from jetset.plot_sedfit import PlotSED
from jetset.test_data_helper import  test_SEDs

sed_data=ObsData.load('Mrk_401.pkl')

ms=McmcSampler.load('mcmc_sampler.pkl')

ms.model.name

'Only-Synch-best-fit-minuit'

p=ms.plot_model(sed_data=sed_data,fit_range=[1E15, 1E20],size=50)
p.setlim(y_min=1E-14,x_min=1E6,x_max=2E20)
../../../_images/Jet_example_only_synchrotron_48_0.png 
p=ms.plot_model(sed_data=sed_data,fit_range=[1E15, 1E20],size=100,quantiles=[0.05,0.95])
p.setlim(y_min=1E-14,x_min=1E6,x_max=2E20)
../../../_images/Jet_example_only_synchrotron_49_0.png 
f=ms.plot_par('beam_obj',log_plot=False)
../../../_images/Jet_example_only_synchrotron_50_0.png 
f=ms.plot_par('B',log_plot=True)
../../../_images/Jet_example_only_synchrotron_51_0.png 
f=mcmc.plot_chain('s',log_plot=False)
../../../_images/Jet_example_only_synchrotron_52_0.png 
f=mcmc.corner_plot()
../../../_images/Jet_example_only_synchrotron_53_0.png 
mcmc.samples.shape

(5120, 5)

mcmc.calls_tot

6400

5120/6400

0.8