Workflow of MUTIS
Here we illustrate the workflow of MUTIS (class version) to analyze correlation between two signals. The general workflow is:
We load the data in a suitable way.
We generate the synthetic light curves, which are used to calculate the confidence levels of correlations.
This is an automatic process if the method used for generating them is just sampling (lc_gen_samp), or flux randomization with the same psd (lc_gen_psd). If we simulate them as an stochastic process (lc_gen_ou) as suggested for example in Tavecchio, Bonnoli and Galanti, 2020 (doi:10.1093/mnras/staa2055) we need to carefully choose the parameters for this stochastic process (currently trying to optimize it).
We generate the correlations (the synthetic ones and the real one) and plot them together.
For this we need to manually choose the time binning.
Juan Escudero
[1]:
# I'm doing a lot of changes to MUTIS while writting this, better reload automatically.
%load_ext autoreload
%autoreload 2
[2]:
import numpy as np
import scipy as sp
import pandas as pd
pd.set_option("display.max_columns", None)
from astropy.time import Time
import matplotlib as mplt
import matplotlib.pyplot as plt
import mutis
Load data
Load it with pandas becase I like the way pandas prints tables. Also useful accesing vars as strings!
[3]:
data = {}
[4]:
data['3mm'] = pd.read_csv('data/mm-I.dat', comment='!')
data['3mm']
[4]:
| jyear | I | dI | |
|---|---|---|---|
| 0 | 2007.336684 | 4.038324 | 0.180 |
| 1 | 2007.384127 | 3.664838 | 0.200 |
| 2 | 2007.603228 | 1.487572 | 0.090 |
| 3 | 2007.665543 | 1.245085 | 0.070 |
| 4 | 2007.714794 | 0.600193 | 0.040 |
| ... | ... | ... | ... |
| 100 | 2020.429704 | 0.910855 | 0.052 |
| 101 | 2020.445397 | 1.940285 | 0.111 |
| 102 | 2020.507706 | 1.482400 | 0.111 |
| 103 | 2020.650145 | 2.236022 | 0.112 |
| 104 | 2020.750106 | 1.926506 | 0.113 |
105 rows × 3 columns
[5]:
data['gamma'] = pd.read_csv('data/gamma-I.dat', comment='!')
data['gamma']
[5]:
| jyear | CFlux | CFluxErr | |
|---|---|---|---|
| 0 | 2008.605065 | 1.894027e-07 | 3.326300e-08 |
| 1 | 2008.624230 | 5.572932e-07 | 4.754800e-08 |
| 2 | 2008.643395 | 5.439230e-07 | 5.684800e-08 |
| 3 | 2008.662560 | 6.763391e-07 | 6.980500e-08 |
| 4 | 2008.681725 | 9.020173e-07 | 7.435100e-08 |
| ... | ... | ... | ... |
| 631 | 2020.763861 | 1.119233e-07 | 5.747500e-08 |
| 632 | 2020.783026 | 9.138770e-08 | 3.188800e-08 |
| 633 | 2020.802191 | 6.064478e-08 | 1.657600e-08 |
| 634 | 2020.821356 | 4.770228e-08 | 2.429000e-08 |
| 635 | 2020.840520 | 1.733534e-07 | 3.886500e-08 |
636 rows × 3 columns
Show data
[6]:
fig, ax = plt.subplots(nrows=2, figsize=(12,4), sharex=True, gridspec_kw={'hspace':0.00})
ax[0].errorbar(data['gamma']['jyear'], data['gamma']['CFlux'], data['gamma']['CFluxErr'], fmt='b.--', label='Gamma-ray flux (0.1 to 200 GeV)', lw=1) # freq?
ax[0].set_ylabel('$\gamma/\mathrm{cm}^2/\mathrm{s}?$')
ax[0].legend()
ax[0].grid()
ax[1].errorbar(data['3mm']['jyear'], data['3mm']['I'], data['3mm']['dI'], fmt='k.--', label='mm', lw=1) #data['3mm'].plot(x='year', y='I')
ax[1].set_ylabel('Jy')
ax[1].legend()
ax[1].grid()
fig.suptitle('Test Source Flux Density')
plt.xlabel('years')
plt.show()
Correlate data
Choose data and methods
Choose the signals to correlate
[7]:
t1, s1, ds1 = data['3mm']['jyear'], data['3mm']['I'], data['3mm']['dI']
t2, s2, ds2 = data['gamma']['jyear'], 1e6*data['gamma']['CFlux'], 1e6*data['gamma']['CFluxErr'] ### ¡Use sensible units!
Create signals with the corresponding data and choose the method for generating the light curves.
[8]:
sig1 = mutis.Signal(t1, s1, ds1, 'lc_gen_psd_nft')
sig2 = mutis.Signal(t2, s2, ds2, 'lc_gen_ou')
Create the correlation with the corresponding signals and choose the method for generating the correltions.
[9]:
correlation = mutis.Correlation(sig1, sig2, 'welsh')
Generate synthetic light curves
HUMAN INTERVENTION
Chose good parameters for time binning!
Good parameters will have a high and similar \(n_i\) for all bins, with bins covering most of the time range.
[10]:
correlation.gen_times(dtmin=0.25)
[11]:
correlation.plot_times()
HUMAN INTERVENTION
If the method for generating the synthetic light curves is lc_gen_ou we need to chose a good set of parameters.
In this case we have set only sig2’s method to be lc_gen_ou because lc_gen_psd is more or less okey for sig1.
The Signal class includes methods to deal with this OU parameter extraction. This could be more or less an authomathic process, but for some (not a negible number of) cases degenerancy of these OU parameter makes human intervention necessary. We need: 1. To plot the histogram and check if the supposed probability density function (pdf) fits nicely (several ways: curve fitting or MLE). 2. To check if the stimates are reasonable and produce nice light curves. 3. Actually generate the light curves.
Plot the histogram and get extracted parameters
[12]:
%%time
fits = sig2.OU_fit()
/home/docs/checkouts/readthedocs.org/user_builds/mutis/envs/latest/lib/python3.7/site-packages/mutis/signal.py:252: RuntimeWarning: invalid value encountered in double_scalars
/ xx ** (ll + 2)
CPU times: user 16.8 s, sys: 36.9 ms, total: 16.9 s
Wall time: 16.8 s
/home/docs/checkouts/readthedocs.org/user_builds/mutis/envs/latest/lib/python3.7/site-packages/mutis/signal.py:252: RuntimeWarning: divide by zero encountered in true_divide
/ xx ** (ll + 2)
/home/docs/checkouts/readthedocs.org/user_builds/mutis/envs/latest/lib/python3.7/site-packages/mutis/signal.py:252: RuntimeWarning: invalid value encountered in true_divide
/ xx ** (ll + 2)
[13]:
fits
[13]:
{'sigma_est': 5.94771645934871,
'curve_fit': (array([1.00758681, 0.12887555]),
array([0.17348078, 0.01219195])),
'MLE_fit': (0.5523802491183334, 0.14125117915448523),
'th_est1': 9.77031709753,
'th_est2': 17.821858473859322}
Check if this parameters generate nice light curves.
[14]:
#%matplotlib inline
sig2.check_gen(fgen='lc_gen_ou', fgen_params=dict(theta=17.8, mu=0.14, sigma=5.9))
Original vs Synthethic:
mean: 0.11241061716509433 / 0.11872079739831864
std: 0.14036970709552712 / 0.1692722331611489
[14]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
Set the parameters
[15]:
sig2.OU_theta, sig2.OU_mu, sig2.OU_sigma = 17.8, 0.13, 5.9
Actually generate the synthetic light curves
[16]:
%%time
correlation.gen_synth(400)
CPU times: user 1.58 s, sys: 0 ns, total: 1.58 s
Wall time: 1.58 s
Generate and plot correlations
[17]:
%%time
correlation.gen_corr();
CPU times: user 7.91 s, sys: 0 ns, total: 7.91 s
Wall time: 7.91 s
[18]:
correlation.plot_corr()
[18]:
<AxesSubplot:>
Find the peaks
[19]:
pd.DataFrame(correlation.peak_find())
[19]:
| x | s | y | signf1s | signif_percent | |
|---|---|---|---|---|---|
| 0 | 0.126204 | 0.128975 | 0.635940 | 7.398281 | 100.00 |
| 1 | 1.388245 | 0.128975 | 0.073888 | 1.025973 | 85.25 |
| 2 | 2.397877 | 0.128975 | -0.251714 | -3.049960 | 0.25 |
| 3 | -2.650286 | 0.128975 | -0.109054 | -1.255481 | 8.75 |
| 4 | -3.155102 | 0.128975 | -0.122123 | -1.358888 | 6.75 |
| 5 | -4.164735 | 0.128975 | -0.288604 | -3.042922 | 0.00 |
| 6 | -5.174367 | 0.128975 | -0.223148 | -2.639614 | 1.75 |
| 7 | 5.174367 | 0.128975 | -0.003976 | -0.037662 | 43.75 |
| 8 | 6.184000 | 0.128975 | 0.117359 | 1.147013 | 86.00 |
| 9 | -6.941224 | 0.128975 | 0.653198 | 6.733913 | 100.00 |
| 10 | 7.193632 | 0.128975 | 0.563514 | 4.479770 | 100.00 |
| 11 | -8.455673 | 0.128975 | 0.315707 | 3.368703 | 99.75 |
| 12 | 9.297034 | 0.128975 | 0.054819 | 0.535690 | 71.75 |
[ ]: