カテゴリー: 未分類

SSC model optimization to blazar SED data using MCMC.

Download files: yysscfit.tgz

Readme.txt

0. Environment 

The following scripts were checked in:
Python 3.7
numpy 1.15.1
scipy 1.1.0
Cython 0.29.14

1. Files
 
yysscfit.py : Python script for the MCMC optimization of the SSC model
SSC_calc_cy.pyx : Cython script for the SSC calculation
setup.py : Python script to build "SSC_calc_cy"
sample.py : Python script to use "yysscfit"
demo_sed.dat : SED data
 
2. Build
 
Only "SSC_calc_cy" needs to be built.
 
$ cython SSC_calc_cy.pyx
$ python3 setup.py build_ext --inplace
 
"SSC_calc_cy.*.so" will be generated in the same directory.
It will be imported in "yysscfit.py".
Depending on versions, some warnings may be generated when "setup.py" is running.
 
3. Tutorial
 
$ python3 sample.py > mcmc_sample.dat
 
If "sample.py" is
 
yy.yysscfit(infile="MJD55252data_logpar.dat", usecols=[0,2,3], n_iteration=10, z=0.031, random_seed=32)

then, "mcmc_sample.dat" should be like:
 
 1 0.001297101984483134 18.634687576670544 957585.1905420072 2.6521412512098218e+44 27690.051953020986 2.2695156866513972 4.855982220271184 100.0 100000000.0 1.0 -2714.8144130083415
 2 0.001245117674215914 17.081115185147606 1278381.5943068536
 1.3612744844185282e+44 40566.313235642614 2.295914566732403
 4.370741791615346 100.0 100000000.0 1.0 -2713.2264705466346
 
 ...

 9 0.0009052483786278515 18.80058920839121 1037602.849494839 5.120608165846359e+43 454918.0986365516 2.5563176245183143 4.5988153641060245 100.0 100000000.0 0.4444444444444444 -1503.177949138101
 10 0.0009052483786278515 18.80058920839121 1037602.849494839 5.120608165846359e+43 454918.0986365516 2.5563176245183143 4.5988153641060245 100.0 100000000.0 0.4 -1503.177949138101
 Time: 3.821684 sec (10 steps)

The columnes are, 1) Iteration, 2) magnetic field (B), 3) Doppler factor (delta_D), 4) Time scale (T), 5) Normalization of the electron energy distribution (K), 6) gamma_b, 7) p_1, 8) p_2, 9) gamma_min, 10) gamma_max, 11) acceptance ratio, and 12) log posterior.
 
This examples only stop MCMC after 10 iterations.
You can make the trace plots and posterior distribution of the SSC parameters from "mcmc_sample.dat" with more iterations (>10^4).
 
4. Descriptions
 
4-1. yysscfit.py
 
Python script for generating the MCMC samples of the SSC parameters.
The script was originally developed by Yurika Yamada, and modified by Makoto Uemura.
 
yysscfit(infile,usecols,n_iteration,z,
              prior_var=None, prior_mu=None, prior_sig=None,
              fixed_var=None, fixed_value=None,
              eed_single_powerlaw=False,
              prm_init= [np.log10(0.001),np.log10(20),np.log10(1361664),np.log10(1e43),np.log10(10**4.5),2.0,4.0,np.log10(10**2),np.log10(10**8)],
              cov_init=[0.04,0.008,0.02,0.6,0.15,0.2,0.2,4.0,7.0],
              nu_sync_ssc=1e19, max_p1=20.0, logr1_fix=2.0,
              logr2_fix=8.0, random_seed=None)
 
Parameters
----------
infile : str
     Name of the input data file including nu, nu_f_nu, and errors.
usecols : list
     List of the column indexes (zero-based) for nu, nu_f_nu, and errors. 
     If the 1st, 2nd, and 3rd columns in "infile" corresponds to nu, nu_f_nu, and errors,    
     then usecols=[0,1,2].
n_iteration : int
     The number of iterations of MCMC. 
     Small number (10--100) should be set for a test.
     A full estimaten requires typically, n_iteration=100000.
z : float
     Redshift
prior_var : str, optional
     Variable of which Gaussian prior is set. 
     It should be either 'logB', 'logD', 'logT', or 'logK'.
     If it is set, you have to give values of prior_mu and prior_sig.
     (Default: None)
prior_mu : float, optional
     The mean of the Gaussian prior in log scale.
prior_sig : float, optional
     The standard deviation of the Gaussian prior in log scale.
fixed_var : str, optional
     Variable fixed. 
     It should be either 'logB', 'logD', 'logT', or 'logK'.
     If it is set, you have to give values of fixed_value.
     (Default: None)
fixed_value : float, optional
     The value for the fixed variable.
eed_single_powerlaw : bool, optional
     If this flag is True, then the model of the electron energy 
     distribution is a single power-law form with a large p1 value 
     of "max_p1". Otherwise, it is a broken power-law form. 
     (Default: False)
prm_init : list, optional
     Initial SSC parameters. The list should include the 9 variables 
     like, [logB, logD, logT, logK, logLC, p0, p1, logr1, logr2]. 
cov_init : list, optional
     Initial values of variances for the covariant matrix of 
     the proposal distribution for MCMC. As in "prm_init", it should 
     be a list having 9 values.
nu_sync_ssc : float, optional
     Frequency below which the synchrotron emission is dominant.
     This code assumes it is dominant below nu_sync_ssc and 
     the SSC is dominant above nu_sync_ssc. It calculates 
     not both components, but either dominant component at each 
     frequency.  (Dafault: 1e19)
max_p1 : float, optional
     Maximum value of p1.  (Default: 20.0)
logr1_fix : float, optional
     Fixed value of log gamma_min. (Default: 2.0)
logr2_fix : float, optional
     Fixed value of log gamma_max. (Default: 8.0)
random_seed : int, optional
     Ramdom seed number. If a value is set, reproductive MCMC samples 
     are obtained.  (Default: None)
 
Returns
-------
 This function only gives outputs to stdout. 
 It outputs a set of SSC parameters for each MCMC iteration.
 It has 12 columns:
 Iteration number, B, D, T, K, LC, p0, p1, r1, r2, acceptance ratio, log posterior
 
 Using the output of "yysscfit", you can make the trace plots and
 posterior distributions of the model parameters.
 "yysscfit" does not provide plotting tools for them.
 
4-2. SSC_calc_cy.pyx & setup.py
 
 Cython script to calcualte synchrotron and SSC flux.
 It was originally developed by Ryosuke Itoh as a python script,
 and modified as a cython script by Makoto Uemura.
 
 If "SSC_calc_cy.pyx" is modified, you need to rebuild it
 by using "setup.py".
 
4-3 sample.py
 
 Python script to use "yysscfit".
 It contains several examples with short descriptions.