.. _model_fitting_2: Model fitting 2: SSC + galaxy template ====================================== .. code:: ipython3 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 .. code:: ipython3 print(jetset.__version__) .. parsed-literal:: 1.3.0rc7 .. code:: ipython3 test_SEDs .. parsed-literal:: ['/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/test_data/SEDs_data/SED_3C345.ecsv', '/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk421_EBL_DEABS.ecsv', '/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv', '/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_DEABS.ecsv'] Loading data ------------ see the :ref:`data_format` user guide for further information about loading data .. code:: ipython3 data=Data.from_file(test_SEDs[3]) .. code:: ipython3 %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() .. parsed-literal:: ================================================================================ *** binning data *** ---> N bins= 90 ---> bin_widht= 0.2 msk [False True False True True True True True False False False True False False False False False False False False False False False False False True True True True True True False False False False False False False True True True True False True True True True True True False True False False False False False False False False False False False False False False False True False True False True False True False True False True False False False False False True True True True True True True False] ================================================================================ .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_8_1.png .. code:: ipython3 sed_data.save('Mrk_501.pkl') Phenomenological model constraining ----------------------------------- see the :ref:`phenom_constr` user guide for further information about loading data Spectral indices ~~~~~~~~~~~~~~~~ .. code:: ipython3 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) .. parsed-literal:: ================================================================================ *** evaluating spectral indices for data *** ================================================================================ .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_13_1.png Sed shaper ~~~~~~~~~~ .. code:: ipython3 mm,best_fit=my_shape.sync_fit(check_host_gal_template=True, Ep_start=None, minimizer='lsb', silent=True, fit_range=[10. , 21.]) .. parsed-literal:: ================================================================================ *** Log-Polynomial fitting of the synchrotron component *** ---> first blind fit run, fit range: [10.0, 21.0] ---> class: HSP ---> class: HSP .. raw:: html Table length=6
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-6.522794e-02-6.522794e-025.892905e-03---4.913172e-02-1.000000e+010.000000e+00False
LogCubicc-1.908748e-03-1.908748e-038.488797e-04--5.440153e-03-1.000000e+011.000000e+01False
LogCubicEp1.704833e+011.704833e+016.858392e-02--1.593204e+010.000000e+003.000000e+01False
LogCubicSp-1.030052e+01-1.030052e+011.424853e-02---1.022242e+01-3.000000e+010.000000e+00False
host_galaxynuFnu_p_host-1.008538e+01-1.008538e+012.900917e-02---1.022242e+01-1.222242e+01-8.222416e+00False
host_galaxynu_scale1.934519e-021.934519e-021.919833e-03--0.000000e+00-5.000000e-015.000000e-01False
.. parsed-literal:: ---> sync nu_p=+1.704833e+01 (err=+6.858392e-02) nuFnu_p=-1.030052e+01 (err=+1.424853e-02) curv.=-6.522794e-02 (err=+5.892905e-03) ================================================================================ .. code:: ipython3 my_shape.IC_fit(fit_range=[23., 29.],minimizer='minuit',silent=True) p=my_shape.plot_shape_fit() p.setlim(y_min=1E-15) .. parsed-literal:: ================================================================================ *** Log-Polynomial fitting of the IC component *** ---> fit range: [23.0, 29.0] ---> LogCubic fit ====> simplex ====> migrad ====> simplex ====> migrad ====> simplex ====> migrad .. raw:: html Table length=4
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-1.323324e-01-1.323324e-013.154579e-02---1.000000e+00-1.000000e+010.000000e+00False
LogCubicc-3.389147e-02-3.389147e-022.022657e-02---1.000000e+00-1.000000e+011.000000e+01False
LogCubicEp2.550950e+012.550950e+012.190510e-01--2.555059e+010.000000e+003.000000e+01False
LogCubicSp-1.057957e+01-1.057957e+014.198643e-02---1.000000e+01-3.000000e+010.000000e+00False
.. parsed-literal:: ---> IC nu_p=+2.550950e+01 (err=+2.190510e-01) nuFnu_p=-1.057957e+01 (err=+4.198643e-02) curv.=-1.323324e-01 (err=+3.154579e-02) ================================================================================ .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_16_3.png Model constraining ~~~~~~~~~~~~~~~~~~ In this step we are not fitting the model, we are just obtaining the phenomenological ``pre_fit`` model, that will be fitted in using minuit ore least-square bound, as shown below .. code:: ipython3 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=1E11, SEDShape=my_shape) prefit_jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=False,silent=True) prefit_jet.save_model('prefit_jet_gal_templ.pkl') .. parsed-literal:: ================================================================================ *** constrains parameters from observable *** ===> setting C threads to 12 .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm1.141280e+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+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*1.487509e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*2.310708e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm32.312656e+010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*1.107634e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.248426e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature3.261397e-01-1.500000e+011.500000e+01FalseFalse
.. parsed-literal:: ================================================================================ .. code:: ipython3 pl=prefit_jet.plot_model(sed_data=sed_data) pl.add_residual_plot(prefit_jet,sed_data) pl.setlim(y_min=1E-14,x_min=1E7,x_max=1E29) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_20_0.png Model fitting procedure ----------------------- .. note:: Please, read the introduction and the caveat :ref:`for the frequentist model fitting ` to understand the frequentist fitting workflow see the :ref:`composite_models` user guide for further information about the implementation of :class:`.FitModel`, in particular for parameter setting Model fitting with LSB ~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 from jetset.model_manager import FitModel from jetset.jet_model import Jet jet=Jet.load_model('prefit_jet_gal_templ.pkl') jet.set_gamma_grid_size(200) .. parsed-literal:: ===> setting C threads to 12 .. code:: ipython3 fit_model=FitModel( jet=jet, name='SSC-best-fit-lsb',template=my_shape.host_gal) fit_model.show_model() .. parsed-literal:: -------------------------------------------------------------------------------- Composite model description -------------------------------------------------------------------------------- name: SSC-best-fit-lsb type: composite_model components models: -model name: jet_leptonic model type: jet -model name: host_galaxy model type: template -------------------------------------------------------------------------------- individual component description -------------------------------------------------------------------------------- model description: -------------------------------------------------------------------------------- type: Jet name: jet_leptonic geometry: spherical electrons distribution: type: lppl gamma energy grid size: 201 gmin grid : 1.487509e+02 gmax grid : 2.310708e+06 normalization: True log-values: False ratio of cold protons to relativistic electrons: 1.000000e+00 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:Sum, hidden: False name:Sync, state: self-abs name:Sync, hidden: False name:SSC, state: on name:SSC, hidden: False 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 -------------------------------------------------------------------------------- .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*1.487509e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*2.310708e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm32.312656e+010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*1.107634e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.248426e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature3.261397e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm1.141280e+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+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseFalse
.. parsed-literal:: -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- model description -------------------------------------------------------------------------------- name: host_galaxy type: template -------------------------------------------------------------------------------- .. raw:: html Table length=2
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
host_galaxynuFnu_p_hostnuFnu-scaleerg / (s cm2)-1.008538e+01-2.000000e+012.000000e+01FalseFalse
host_galaxynu_scalenu-scaleHz1.934519e-02-2.000000e+002.000000e+00FalseTrue
.. parsed-literal:: -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- .. note:: Since the `jet_leptonic to model` has to be summed to the `host_galaxy` model, we do not need to define the functional form for the composite model, because the default compostion is the sum of all the components (see the :ref:`composite_models` user guide for further information about the new implementation of `FitModel`, in particular for parameter setting). Anyhow, we show here the definition of the model composition for purpose of clarity .. code:: ipython3 fit_model.composite_expr='jet_leptonic + host_galaxy' .. code:: ipython3 fit_model.freeze('jet_leptonic','z_cosm') fit_model.freeze('jet_leptonic','R_H') fit_model.jet_leptonic.parameters.N.fit_range=[1E-5, 1E5] fit_model.jet_leptonic.parameters.B.fit_range=[1E-3, 1] fit_model.jet_leptonic.parameters.beam_obj.fit_range=[5., 50.] fit_model.jet_leptonic.parameters.R.fit_range=[1E15,1E17] fit_model.jet_leptonic.parameters.gmax.fit_range=[1E4,1E8] fit_model.host_galaxy.parameters.nuFnu_p_host.frozen=False fit_model.host_galaxy.parameters.nu_scale.frozen=True .. code:: ipython3 from jetset.minimizer import fit_SED,ModelMinimizer model_minimizer_lsb=ModelMinimizer('lsb') best_fit_lsb=model_minimizer_lsb.fit(fit_model,sed_data,1E11,1E29,fitname='SSC-best-fit-lsb',repeat=1) .. parsed-literal:: filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 ================================================================================ *** start fit process *** ----- .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: - best chisq=4.16923e+01 ------------------------------------------------------------------------- Fit report Model: SSC-best-fit-lsb .. raw:: html Table length=14
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*8.580927e+011.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*2.310708e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm35.244887e+010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*1.069947e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.220278e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature2.813492e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm1.141280e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss3.716188e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.575826e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseTrue
host_galaxynuFnu_p_hostnuFnu-scaleerg / (s cm2)-1.005821e+01-2.000000e+012.000000e+01FalseFalse
host_galaxynu_scalenu-scaleHz1.934519e-02-2.000000e+002.000000e+00FalseTrue
.. parsed-literal:: converged=True calls=429 mesg= .. parsed-literal:: '`ftol` termination condition is satisfied.' .. parsed-literal:: dof=21 chisq=41.692309, chisq/red=1.985348 null hypothesis sig=0.004599 best fit pars .. raw:: html Table length=14
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
jet_leptonicgmin8.580927e+018.580927e+012.245405e+02--1.487509e+021.000000e+001.000000e+09False
jet_leptonicgmax2.310708e+062.310708e+063.063708e+06--2.310708e+061.000000e+041.000000e+08False
jet_leptonicN5.244887e+015.244887e+013.642249e+02--2.312656e+011.000000e-051.000000e+05False
jet_leptonicgamma0_log_parab1.069947e+041.069947e+041.593138e+04--1.107634e+041.000000e+001.000000e+09False
jet_leptonics2.220278e+002.220278e+006.119317e-02--2.248426e+00-1.000000e+011.000000e+01False
jet_leptonicr2.813492e-012.813492e-016.199616e-02--3.261397e-01-1.500000e+011.500000e+01False
jet_leptonicR1.141280e+161.141280e+164.323108e+16--1.141280e+161.000000e+151.000000e+17False
jet_leptonicR_H1.000000e+17------1.000000e+170.000000e+00--True
jet_leptonicB3.716188e-023.716188e-024.850783e-02--5.050000e-021.000000e-031.000000e+00False
jet_leptonicNH_cold_to_rel_e1.000000e+00------1.000000e+000.000000e+00--True
jet_leptonicbeam_obj2.575826e+012.575826e+013.748732e+01--2.500000e+015.000000e+005.000000e+01False
jet_leptonicz_cosm3.360000e-02------3.360000e-020.000000e+00--True
host_galaxynuFnu_p_host-1.005821e+01-1.005821e+013.476096e-02---1.008538e+01-1.222242e+01-8.222416e+00False
host_galaxynu_scale1.934519e-02------1.934519e-02-5.000000e-015.000000e-01True
.. parsed-literal:: ------------------------------------------------------------------------- ================================================================================ .. code:: ipython3 best_fit_lsb.save_report('SSC-best-fit-lsb.pkl') model_minimizer_lsb.save_model('model_minimizer_lsb.pkl') fit_model.save_model('fit_model_lsb.pkl') .. code:: ipython3 %matplotlib inline fit_model.set_nu_grid(1E6,1E30,200) fit_model.eval() p2=fit_model.plot_model(sed_data=sed_data) p2.setlim(y_min=1E-14,x_min=1E6,x_max=2E28) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_31_0.png Model fitting with Minuit ~~~~~~~~~~~~~~~~~~~~~~~~~ To run the ``minuit`` minimizer we will use the best-fit results from ``lsb`` to set the boundaries for our parameters. .. code:: ipython3 fit_model.freeze('jet_leptonic','z_cosm') fit_model.freeze('jet_leptonic','R_H') fit_model.jet_leptonic.parameters.beam_obj.fit_range=[5., 50.] fit_model.jet_leptonic.parameters.R.fit_range=[10**15.5,10**17.5] fit_model.host_galaxy.parameters.nuFnu_p_host.frozen=False fit_model.host_galaxy.parameters.nu_scale.frozen=True fit_model.jet_leptonic.parameters.gmin.fit_range=[10,1000] fit_model.jet_leptonic.parameters.gmax.fit_range=[5E5,1E8] fit_model.jet_leptonic.parameters.gamma0_log_parab.fit_range=[1E3,5E5] .. code:: ipython3 model_minimizer_minuit=ModelMinimizer('minuit') best_fit_minuit=model_minimizer_minuit.fit(fit_model,sed_data,1E11,1E29,fitname='SSC-best-fit-minuit',repeat=3) .. parsed-literal:: filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 ================================================================================ *** start fit process *** ----- fit run: 0 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> simplex ====> migrad - best chisq=1.04513e+01 fit run: 1 - old chisq=1.04513e+01 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> simplex ====> migrad - best chisq=1.04509e+01 fit run: 2 - old chisq=1.04509e+01 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> simplex ====> migrad - best chisq=1.04509e+01 ------------------------------------------------------------------------- Fit report Model: SSC-best-fit-minuit .. raw:: html Table length=14
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*5.455560e+011.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*2.105148e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm34.504764e+010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*5.386110e+031.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.168634e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature2.323139e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm1.298307e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.206137e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming4.763830e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseTrue
host_galaxynuFnu_p_hostnuFnu-scaleerg / (s cm2)-1.008451e+01-2.000000e+012.000000e+01FalseFalse
host_galaxynu_scalenu-scaleHz1.934519e-02-2.000000e+002.000000e+00FalseTrue
.. parsed-literal:: converged=True calls=1019 mesg= .. raw:: html
Migrad
FCN = 10.45 Nfcn = 1019
EDM = 1.19e-05 (Goal: 0.0002) time = 27.0 sec
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 par_0 54.56 0.27 10 1E+03
1 par_1 2.11e6 0.17e6 5E+05 1E+08
2 par_2 45.0 0.4 1E-05 1E+05
3 par_3 5.39e3 0.08e3 1E+03 5E+05
4 par_4 2.169 0.004 -10 10
5 par_5 232.3e-3 0.8e-3 -15 15
6 par_6 12.983e15 0.033e15 3.16E+15 3.16E+17
7 par_7 12.06e-3 0.09e-3 0.001 1
8 par_8 47.64 0.17 5 50
9 par_9 -10.085 0.020 -12.2 -8.22
.. parsed-literal:: dof=21 chisq=10.450851, chisq/red=0.497660 null hypothesis sig=0.972444 best fit pars .. raw:: html Table length=14
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
jet_leptonicgmin5.455560e+015.455560e+012.726560e-01--8.580927e+011.000000e+011.000000e+03False
jet_leptonicgmax2.105148e+062.105148e+061.675281e+05--2.310708e+065.000000e+051.000000e+08False
jet_leptonicN4.504764e+014.504764e+014.131515e-01--5.244887e+011.000000e-051.000000e+05False
jet_leptonicgamma0_log_parab5.386110e+035.386110e+037.892077e+01--1.069947e+041.000000e+035.000000e+05False
jet_leptonics2.168634e+002.168634e+003.851834e-03--2.220278e+00-1.000000e+011.000000e+01False
jet_leptonicr2.323139e-012.323139e-017.857833e-04--2.813492e-01-1.500000e+011.500000e+01False
jet_leptonicR1.298307e+161.298307e+163.259461e+13--1.141280e+163.162278e+153.162278e+17False
jet_leptonicR_H1.000000e+17------1.000000e+170.000000e+00--True
jet_leptonicB1.206137e-021.206137e-028.551352e-05--3.716188e-021.000000e-031.000000e+00False
jet_leptonicNH_cold_to_rel_e1.000000e+00------1.000000e+000.000000e+00--True
jet_leptonicbeam_obj4.763830e+014.763830e+011.659000e-01--2.575826e+015.000000e+005.000000e+01False
jet_leptonicz_cosm3.360000e-02------3.360000e-020.000000e+00--True
host_galaxynuFnu_p_host-1.008451e+01-1.008451e+011.998494e-02---1.005821e+01-1.222242e+01-8.222416e+00False
host_galaxynu_scale1.934519e-02------1.934519e-02-5.000000e-015.000000e-01True
.. parsed-literal:: ------------------------------------------------------------------------- ================================================================================ for further information regardin minuit please refer to https://iminuit.readthedocs.io/en/stable/ .. code:: ipython3 %matplotlib inline fit_model.set_nu_grid(1E6,1E30,200) fit_model.eval() p2=fit_model.plot_model(sed_data=sed_data) p2.setlim(y_min=1E-14,x_min=1E6,x_max=2E28) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_37_0.png .. code:: ipython3 best_fit_minuit.save_report('SSC-best-fit-minuit.pkl') model_minimizer_minuit.save_model('model_minimizer_minuit.pkl') fit_model.save_model('fit_model_minuit.pkl') Model fitting with a bkn pl ~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 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='bkn', t_var_sec=3*86400, nu_cut_IR=1E11, SEDShape=my_shape) prefit_jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=False,silent=True) prefit_jet.save_model('prefit_jet_bkn_gal_templ.pkl') .. parsed-literal:: ================================================================================ *** constrains parameters from observable *** ===> setting C threads to 12 .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm1.247372e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss2.847716e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*1.980875e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*3.077106e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.406682e+010.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*2.094216e+051.000000e+001.000000e+09FalseFalse
jet_leptonicpLE_spectral_slope2.248426e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicp_1HE_spectral_slope3.500000e+00-1.000000e+011.000000e+01FalseFalse
.. parsed-literal:: ================================================================================ .. code:: ipython3 pl=prefit_jet.plot_model(sed_data=sed_data) pl.add_residual_plot(prefit_jet,sed_data) pl.setlim(y_min=1E-14,x_min=1E7,x_max=1E29) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_41_0.png .. code:: ipython3 jet_minuit_bkn=Jet.load_model('prefit_jet_bkn_gal_templ.pkl') jet_minuit_bkn.set_gamma_grid_size(200) fit_model_bkn=FitModel( jet=jet_minuit_bkn, name='SSC-best-fit-bkn-lsb',template=my_shape.host_gal) fit_model_bkn.freeze('jet_leptonic','z_cosm') fit_model_bkn.freeze('jet_leptonic','R_H') fit_model_bkn.jet_leptonic.parameters.beam_obj.fit_range=[5,50] fit_model_bkn.jet_leptonic.parameters.R.fit_range=[10**15.5,10**17.5] fit_model_bkn.jet_leptonic.parameters.gmax.fit_range=[1E4,1E8] fit_model_bkn.host_galaxy.parameters.nuFnu_p_host.frozen=False fit_model_bkn.host_galaxy.parameters.nu_scale.frozen=True .. parsed-literal:: ===> setting C threads to 12 .. code:: ipython3 model_minimizer_lsb_bkn=ModelMinimizer('lsb') best_fit_lsb_bkn=model_minimizer_lsb_bkn.fit(fit_model_bkn,sed_data,1E11,1E29,fitname='SSC-best-fit-lsb') .. parsed-literal:: filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 ================================================================================ *** start fit process *** ----- .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: - best chisq=2.49214e+01 ------------------------------------------------------------------------- Fit report Model: SSC-best-fit-lsb .. raw:: html Table length=14
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*1.485242e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*3.071765e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.629672e+010.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*8.193745e+041.000000e+001.000000e+09FalseFalse
jet_leptonicpLE_spectral_slope2.262508e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicp_1HE_spectral_slope3.169331e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicRregion_sizecm1.247372e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.430585e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming4.463828e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseTrue
host_galaxynuFnu_p_hostnuFnu-scaleerg / (s cm2)-1.006917e+01-2.000000e+012.000000e+01FalseFalse
host_galaxynu_scalenu-scaleHz1.934519e-02-2.000000e+002.000000e+00FalseTrue
.. parsed-literal:: converged=True calls=482 mesg= .. parsed-literal:: '`ftol` termination condition is satisfied.' .. parsed-literal:: dof=21 chisq=24.921362, chisq/red=1.186732 null hypothesis sig=0.250586 best fit pars .. raw:: html Table length=14
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
jet_leptonicgmin1.485242e+021.485242e+021.195436e+02--1.980875e+021.000000e+001.000000e+09False
jet_leptonicgmax3.071765e+063.071765e+061.921033e+06--3.077106e+061.000000e+041.000000e+08False
jet_leptonicN1.629672e+011.629672e+012.985576e+01--1.406682e+010.000000e+00--False
jet_leptonicgamma_break8.193745e+048.193745e+045.261749e+04--2.094216e+051.000000e+001.000000e+09False
jet_leptonicp2.262508e+002.262508e+003.933694e-02--2.248426e+00-1.000000e+011.000000e+01False
jet_leptonicp_13.169331e+003.169331e+004.793906e-02--3.500000e+00-1.000000e+011.000000e+01False
jet_leptonicR1.247372e+161.247372e+162.236984e+16--1.247372e+163.162278e+153.162278e+17False
jet_leptonicR_H1.000000e+17------1.000000e+170.000000e+00--True
jet_leptonicB1.430585e-021.430585e-021.084428e-02--2.847716e-020.000000e+00--False
jet_leptonicNH_cold_to_rel_e1.000000e+00------1.000000e+000.000000e+00--True
jet_leptonicbeam_obj4.463828e+014.463828e+012.806607e+01--2.500000e+015.000000e+005.000000e+01False
jet_leptonicz_cosm3.360000e-02------3.360000e-020.000000e+00--True
host_galaxynuFnu_p_host-1.006917e+01-1.006917e+012.666541e-02---1.008451e+01-1.222242e+01-8.222416e+00False
host_galaxynu_scale1.934519e-02------1.934519e-02-5.000000e-015.000000e-01True
.. parsed-literal:: ------------------------------------------------------------------------- ================================================================================ .. code:: ipython3 %matplotlib inline fit_model_bkn.set_nu_grid(1E6,1E30,200) fit_model_bkn.eval() p2=fit_model_bkn.plot_model(sed_data=sed_data) p2.setlim(y_min=1E-14,x_min=1E6,x_max=2E28) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_44_0.png .. code:: ipython3 fit_model_bkn.jet_leptonic.parameters.beam_obj.fit_range=[5,50] fit_model_bkn.jet_leptonic.parameters.R.fit_range=[10**15.5,10**17.5] fit_model_bkn.host_galaxy.parameters.nuFnu_p_host.frozen=False fit_model_bkn.host_galaxy.parameters.nu_scale.frozen=True fit_model_bkn.jet_leptonic.parameters.gmin.fit_range=[10,1000] fit_model_bkn.jet_leptonic.parameters.gmax.fit_range=[5E5,1E8] fit_model_bkn.jet_leptonic.parameters.gamma_break.fit_range=[1E3,1E6] fit_model_bkn.jet_leptonic.parameters.p.fit_range=[1,3] fit_model_bkn.jet_leptonic.parameters.p_1.fit_range=[2,5] model_minimizer_minuit_bkn=ModelMinimizer('minuit') best_fit_minuit_bkn=model_minimizer_minuit.fit(fit_model_bkn,sed_data,1E11,1E29,fitname='SSC-best-fit-minuit-bkn',repeat=3) .. parsed-literal:: filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 ================================================================================ *** start fit process *** ----- fit run: 0 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> simplex ====> migrad - best chisq=1.34228e+01 fit run: 1 - old chisq=1.34228e+01 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> simplex ====> migrad - best chisq=1.30362e+01 fit run: 2 - old chisq=1.30362e+01 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> simplex ====> migrad - best chisq=1.30275e+01 ------------------------------------------------------------------------- Fit report Model: SSC-best-fit-minuit-bkn .. raw:: html Table length=14
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*1.345014e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.773731e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.625170e+010.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*5.531173e+041.000000e+001.000000e+09FalseFalse
jet_leptonicpLE_spectral_slope2.248959e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicp_1HE_spectral_slope2.952674e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicRregion_sizecm1.464854e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.349383e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming4.128886e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseTrue
host_galaxynuFnu_p_hostnuFnu-scaleerg / (s cm2)-1.007170e+01-2.000000e+012.000000e+01FalseFalse
host_galaxynu_scalenu-scaleHz1.934519e-02-2.000000e+002.000000e+00FalseTrue
.. parsed-literal:: converged=True calls=1353 mesg= .. raw:: html
Migrad
FCN = 13.03 Nfcn = 1353
EDM = 5.33e-06 (Goal: 0.0002) time = 27.0 sec
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 par_0 134.501 0.019 10 1E+03
1 par_1 1.7737e6 0.0011e6 5E+05 1E+08
2 par_2 16.2517 0.0009 0
3 par_3 55.312e3 0.032e3 1E+03 1E+06
4 par_4 2.248959 0.000020 1 3
5 par_5 2.95267 0.00015 2 5
6 par_6 14.65e15 0.05e15 3.16E+15 3.16E+17
7 par_7 13.494e-3 0.010e-3 0
8 par_8 41.29 0.11 5 50
9 par_9 -10.072 0.015 -12.2 -8.22
.. parsed-literal:: dof=21 chisq=13.027500, chisq/red=0.620357 null hypothesis sig=0.907658 best fit pars .. raw:: html Table length=14
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
jet_leptonicgmin1.345014e+021.345014e+021.935667e-02--1.485242e+021.000000e+011.000000e+03False
jet_leptonicgmax1.773731e+061.773731e+061.129851e+03--3.071765e+065.000000e+051.000000e+08False
jet_leptonicN1.625170e+011.625170e+018.522217e-04--1.629672e+010.000000e+00--False
jet_leptonicgamma_break5.531173e+045.531173e+043.220423e+01--8.193745e+041.000000e+031.000000e+06False
jet_leptonicp2.248959e+002.248959e+002.045150e-05--2.262508e+001.000000e+003.000000e+00False
jet_leptonicp_12.952674e+002.952674e+001.524221e-04--3.169331e+002.000000e+005.000000e+00False
jet_leptonicR1.464854e+161.464854e+165.023364e+13--1.247372e+163.162278e+153.162278e+17False
jet_leptonicR_H1.000000e+17------1.000000e+170.000000e+00--True
jet_leptonicB1.349383e-021.349383e-021.023312e-05--1.430585e-020.000000e+00--False
jet_leptonicNH_cold_to_rel_e1.000000e+00------1.000000e+000.000000e+00--True
jet_leptonicbeam_obj4.128886e+014.128886e+011.130872e-01--4.463828e+015.000000e+005.000000e+01False
jet_leptonicz_cosm3.360000e-02------3.360000e-020.000000e+00--True
host_galaxynuFnu_p_host-1.007170e+01-1.007170e+011.541061e-02---1.006917e+01-1.222242e+01-8.222416e+00False
host_galaxynu_scale1.934519e-02------1.934519e-02-5.000000e-015.000000e-01True
.. parsed-literal:: ------------------------------------------------------------------------- ================================================================================ .. code:: ipython3 %matplotlib inline fit_model_bkn.set_nu_grid(1E6,1E30,200) fit_model_bkn.eval() p2=fit_model_bkn.plot_model(sed_data=sed_data) p2.setlim(y_min=1E-14,x_min=1E6,x_max=2E28) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_46_0.png .. code:: ipython3 %matplotlib inline from jetset.plot_sedfit import PlotSED fit_model_bkn.set_nu_grid(1E6,1E30,200) fit_model_bkn.eval() fit_model.set_nu_grid(1E6,1E30,200) fit_model.eval() p2=PlotSED() p2.add_data_plot(sed_data,fit_range=[ 1E11, 1E29]) p2.add_model_plot(fit_model,color='black') p2.add_residual_plot(fit_model,sed_data,fit_range=[ 1E11, 1E29],color='black') p2.add_model_plot(fit_model_bkn,color='red') p2.add_residual_plot(fit_model_bkn,sed_data,fit_range=[ 1E11, 1E29],color='red') p2.setlim(y_min=1E-14,x_min=1E6,x_max=2E28) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_47_0.png MCMC sampling ------------- .. note:: Please, read the introduction and the caveat :ref:`for the Bayesian model fitting ` to understand the MCMC sampler workflow. .. code:: ipython3 from jetset.mcmc import McmcSampler from jetset.minimizer import ModelMinimizer .. code:: ipython3 model_minimizer_lsb = ModelMinimizer.load_model('model_minimizer_minuit.pkl') mcmc=McmcSampler(model_minimizer_lsb) .. parsed-literal:: ===> setting C threads to 12 .. code:: ipython3 labels=['N','B','beam_obj','s','gamma0_log_parab'] model_name='jet_leptonic' use_labels_dict={model_name:labels} mcmc.set_labels(use_labels_dict=use_labels_dict) .. code:: ipython3 mcmc.set_bounds(bound=5.0,bound_rel=True) .. parsed-literal:: par: N best fit value: 45.0476406721727 mcmc bounds: [1e-05, 270.2858440330362] par: B best fit value: 0.012061365634053165 mcmc bounds: [0.001, 0.07236819380431898] par: beam_obj best fit value: 47.638299687949946 mcmc bounds: [5.0, 50.0] par: s best fit value: 2.1686339157854624 mcmc bounds: [-8.67453566314185, 10] par: gamma0_log_parab best fit value: 5386.109843065869 mcmc bounds: [1000.0, 32316.659058395213] .. code:: ipython3 mcmc.run_sampler(nwalkers=20, burnin=50,steps=500,progress='notebook') .. parsed-literal:: ===> setting C threads to 12 mcmc run starting .. parsed-literal:: 0%| | 0/500 [00:00 setting C threads to 12 ===> setting C threads to 12 .. code:: ipython3 ms.model.nu_min=1E6 ms.model.jet_leptonic.nu_min=1E6 p=ms.plot_model(sed_data=sed_data,fit_range=[2E11, 2E28],size=100) p.setlim(y_min=1E-14,x_min=1E6,x_max=2E28) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_70_0.png .. code:: ipython3 p=ms.plot_model(sed_data=sed_data,fit_range=[2E11, 2E28],quantiles=[0.05,0.95]) p.setlim(y_min=1E-14,x_min=1E6,x_max=2E28) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_71_0.png .. code:: ipython3 f=ms.plot_par('beam_obj',log_plot=False) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_72_0.png .. code:: ipython3 f=ms.plot_chain(log_plot=False) .. image:: Jet_example_model_fit_wiht_gal_template_files/Jet_example_model_fit_wiht_gal_template_73_0.png