Spectral fit
The pulsar spectral fitting is explained in Swainston et al. 2022 and based on Jankowski et al. 2018. We will summarise how the fitting is done and examples of how to improve your fits.
Fitting algorithm
To account for underestimated uncertainties on outlier points, we modify the regular least-squared quadratic loss function to deviate to linear loss once a certain distance is reached from the model. In this way, outlier data are penalised, and bad data is less likely to skew the model fit. We use the Huber loss function, which is defined as
where \(t\) is a residual and \(k\) is a constant (which we set to 1.345) that defines the point at which outlying points are penalised.
The code will use mingrad function from the iminuit Python package is used to find the minimum of the cost function. This uses a maximum of 10000 calls to converge with the Estimated Distance to Minimum (EDM) criterion. We set the tolerance so that the EDM must be less than \(10^{-8} * \\mathrm{errordef}\).
In the rare cases that migrad does not find a valid fit, we then try the simplex and scan minimisation methods. Simplex does not use derivatives, making it slower but can perform better in some instances. Scan is a brute-force minimisation that uses a grid of possible solutions based on the limits of the input parameters.
The uncertainties were computed using hesse, an error calculator which computes the covariance matrix for the fitted parameters and determines
the \(1\sigma\) uncertainties as the square root of the diagonal elements.
This is all done within the pulsar_spectra.spectral_fit.iminuit_fit_spectral_model() function.
Models
This fit is done for all functions in the models module that are included in pulsar_spectra.models.model_settings().
For example, at the time of writing this documentation, the list of models within model settings includes:
model_dict = {
# Name: [model_function, short_name, start_params, mod_limits]
"simple_power_law" : [
simple_power_law,
"simple pl",
# (a, c)
(a_s, c_s),
[(a_min, a_max), (c_min, c_max)],
simple_power_law_integrate,
],
"broken_power_law" : [
broken_power_law,
"broken pl",
#(vb, a1, a2, c)
(1e9, a_s, a_s, c_s),
[(50e6, 5e9), (a_min, a_max), (a_min, a_max), (c_min, c_max)],
broken_power_law_intergral,
],
"high_frequency_cut_off_power_law" : [
high_frequency_cut_off_power_law,
"pl hard cut-off",
#(vc, a, c)
(vc_s, a_s, c_s),
[vc_both, (a_min, 0.), (c_min, c_max)],
high_frequency_cut_off_power_law_taylor,
],
"low_frequency_turn_over_power_law" : [
low_frequency_turn_over_power_law,
"pl low turn-over",
#(vpeak, a, c, beta)
(vpeak_s, a_s, c_s, beta_s),
[(vpeak_min, vpeak_max), (a_min, 0.), (c_min, c_max) , (beta_min, beta_max)],
low_frequency_turn_over_power_law_taylor,
],
"double_turn_over_spectrum" : [
double_turn_over_spectrum,
"double turn over spectrum",
#(vc, vpeak, a, beta, c)
(vc_s, vpeak_s, a_s, beta_s, c_s),
[(vc_both), (vpeak_min, vpeak_max), (a_min, 0.), (beta_min, beta_max), (c_min, c_max)],
double_turn_over_spectrum_taylor,
],
}
Each item in the dictionary is one of the models that the fitting code will use to fit the pulsar’s spectra. Each item includes a list of the model function, a short name (for plotting), the starting value for each parameter, and the fit limits for each parameter. You can change some of the starting parameters of fit limits if you think it will improve the fit or even comment out a model you do not want to use, like so:
model_dict = {
# Name: [model_function, short_name, start_params, mod_limits]
"simple_power_law" : [
simple_power_law,
"simple pl",
# (a, c)
(a_s, c_s),
[(a_min, a_max), (c_min, c_max)],
simple_power_law_integrate,
],
"broken_power_law" : [
broken_power_law,
"broken pl",
#(vb, a1, a2, c)
(1e9, a_s, a_s, c_s),
[(50e6, 5e9), (a_min, a_max), (a_min, a_max), (c_min, c_max)],
broken_power_law_intergral,
],
"high_frequency_cut_off_power_law" : [
high_frequency_cut_off_power_law,
"pl hard cut-off",
#(vc, a, c)
(vc_s, a_s, c_s),
[vc_both, (a_min, 0.), (c_min, c_max)],
high_frequency_cut_off_power_law_taylor,
],
"low_frequency_turn_over_power_law" : [
low_frequency_turn_over_power_law,
"pl low turn-over",
#(vpeak, a, c, beta)
(vpeak_s, a_s, c_s, beta_s),
[(vpeak_min, vpeak_max), (a_min, 0.), (c_min, c_max) , (beta_min, beta_max)],
low_frequency_turn_over_power_law_taylor,
],
# "double_turn_over_spectrum" : [
# double_turn_over_spectrum,
# "double turn over spectrum",
# #(vc, vpeak, a, beta, c)
# (vc_s, vpeak_s, a_s, beta_s, c_s),
# [(vc_both), (vpeak_min, vpeak_max), (a_min, 0.), (beta_min, beta_max), (c_min, c_max)],
# double_turn_over_spectrum_taylor,
# ],
}
So now, once you reinstall the software, the code will not fit a double turn over spectrum model.
Checking which models you are using
If you are unsure which models or derivations you are using in your fitting, you can use the following function option to print the models info like so:
from pulsar_spectra.models import model_settings
model_settings(print_models=True)
Which will output something like this:
simple_power_law
model_function: simple_power_law
model_function_integrate: simple_power_law_integrate
short_name: simple pl
start_params: (-1.6, 1.0)
mod_limits: [(-8.0, 3.0), (0.0, None)]
broken_power_law
model_function: broken_power_law
model_function_integrate: broken_power_law_intergral
short_name: broken pl
start_params: (1000000000.0, -1.6, -1.6, 1.0)
mod_limits: [(50000000.0, 5000000000.0), (-8.0, 3.0), (-8.0, 3.0), (0.0, None)]
high_frequency_cut_off_power_law
model_function: high_frequency_cut_off_power_law
model_function_integrate: high_frequency_cut_off_power_law_taylor
short_name: pl hard cut-off
start_params: (4000000000.0, -1.6, 1.0)
mod_limits: [None, (-8.0, 0.0), (0.0, None)]
low_frequency_turn_over_power_law
model_function: low_frequency_turn_over_power_law
model_function_integrate: low_frequency_turn_over_power_law_taylor
short_name: pl low turn-over
start_params: (100000000.0, -1.6, 1.0, 1.0)
mod_limits: [(10000000.0, 2000000000.0), (-8.0, 0.0), (0.0, None), (0.1, 2.1)]
double_turn_over_spectrum
model_function: double_turn_over_spectrum
model_function_integrate: double_turn_over_spectrum_taylor
short_name: double turn over spectrum
start_params: (4000000000.0, 100000000.0, -1.6, 1.0, 1.0)
mod_limits: [None, (10000000.0, 2000000000.0), (-8.0, 0.0), (0.1, 2.1), (0.0, None)]
You can find the descriptions of the models in the the models module.
Adding a new model
If you would like to use a new model, you can add a function to the models’ module and set up the defaults for its
initial fit parameters and limits in pulsar_spectra.models.model_settings().
For example, here is the function for the simple power law in the models module:
def simple_power_law(v, a, c, v0):
"""Simple power law:
.. math::
S_v = c \\left( \\frac{v}{v_0} \\right)^a
Parameters
----------
v : `list`
Frequency in Hz.
a : `float`
Spectral Index.
c : `float`
Constant.
v0 : `float`
Reference frequency.
Returns
-------
S_v : `list`
The flux density predicted by the model.
"""
return c*(v/v0)**a
This is the format you must follow to add your model. Frequency must be the first argument, reference frequency must be the last, and we recommend you make a docstring as shown in the above example.
As explained in the previous section, you must add your new model to pulsar_spectra.models.model_settings().
Here are the values for the simple power law:
# fit starting value, min and max
# constant
c_s = 1.
c_min = 0.
c_max = None
# spectral index
a_s = -1.6
a_min = -8.
a_max = 3.
model_dict = {
# Name: [model_function, short_name, start_params, mod_limits]
"simple_power_law" : [
simple_power_law,
"simple pl",
# (a, c)
(a_s, c_s),
[(a_min, a_max), (c_min, c_max)],
],
Because some of the models have common parameters (such as spectral index), some of the fit values have been predefined to be consistent between models.
Make sure you reinstall pulsar_spectra to apply any changes you have made to pulsar_spectra.models.model_settings(), then you will be ready to fit with your new model.
Best fit
The best fit model is determined using the Akaike information criterion (AIC), which measures how much information the model retains about the data without overfitting. It was implemented as
where \(\beta_\mathrm{min}\) is the minimised robust cost function, \(K\) is the number of free parameters, and \(N\) is the number of data points in the fit. The last term is the correction for finite sample sizes, which goes to zero as the sample size gets sufficiently large. The model which results in the lowest AIC is the most likely to be the best fitting model.
All of this is done by calling the pulsar_spectra.spectral_fit.find_best_spectral_fit() function like so:
from pulsar_spectra.catalogue import collect_catalogue_fluxes
from pulsar_spectra.spectral_fit import find_best_spectral_fit
cat_dict = collect_catalogue_fluxes()
pulsar = 'J1453-6413'
freqs, fluxs, flux_errs, refs = cat_dict[pulsar]
best_model_name, iminuit_result, fit_info, p_best, p_category = find_best_spectral_fit(pulsar, freqs, fluxs, flux_errs, refs, plot_best=True)
To confirm that the best model has been found, you can visually inspect the fits of all models using the plot_compare option like so
best_model_name, iminuit_result, fit_info, p_best, p_category = find_best_spectral_fit(pulsar, freqs, fluxs, flux_errs, refs, plot_compare=True)
which will produce
From this plot, it does look like the power-law with a low-frequency turnover is the best model as the code predicted.
If this is not the case and wanted to try and improve the broken power-law fit, for example, you can have more control over the
fit using pulsar_spectra.spectral_fit.iminuit_fit_spectral_model() function like so.
from pulsar_spectra.catalogue import collect_catalogue_fluxes
from pulsar_spectra.spectral_fit import iminuit_fit_spectral_model
cat_list = collect_catalogue_fluxes()
pulsar = 'J1453-6413'
freqs, fluxs, flux_errs, refs = cat_list[pulsar]
# Broken power law function is in the format
# broken_power_law(v, vb, a1, a2, b, v0)
# start params for (v, vb, a1, a2, b)
start_params = (5e8, -1.6, -1.6, 0.1)
# Fit param limits (min, max) or (v, vb, a1, a2, b)
mod_limits = [(None, None), (-10, 10), (-10, 0), (0, None)]
# None means there is no limit
aic, iminuit_result, fit_info = iminuit_fit_spectral_model(
freqs,
fluxs,
flux_errs,
refs,
model_name="broken_power_law",
start_params=start_params,
mod_limits=mod_limits,
plot=True,
save_name="J1453-6413_broken_power_law.png",
)
In this example we are manually handing pulsar_spectra.spectral_fit.iminuit_fit_spectral_model() the default
start_params and mod_limits but you can edit these.