Synthethic generation of light curves
This notebooks illustrates the synthethic generation of light curves with mutis.
There are three methods implemented with mutis, with some variations: - Generating signals by sampling. This signals have same statisical distribution. - Gnerating signals by randomization of the Fourier transform from the PSD (Power Spectral Distribution). This signals have the same variability, and same mean and std. - Generating signals by integration of an stochastic process (Orstein-Uhlenbeck process). Signals have similar shape.
[1]:
# I'm doing a lot of changes to MUTIS while writting this, better reload automatically.
%load_ext autoreload
%autoreload
import numpy as np
import pandas as pd
import matplotlib as mplt
import matplotlib.pyplot as plt
from astropy.time import Time
import mutis
Load test data
[2]:
data = {}
[3]:
data['3mm'] = pd.read_csv('data/mm-I.dat', comment='!')
data['3mm']
[3]:
| 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
[4]:
data['gamma'] = pd.read_csv('data/gamma-I.dat', comment='!')
data['gamma']
[4]:
| 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
3mm
[5]:
import mutis
sig3mm = mutis.Signal(data['3mm']['jyear'], data['3mm']['I'], data['3mm']['dI'])
sig3mm.plot()
[5]:
[<matplotlib.lines.Line2D at 0x7ff878a95f90>]
lc_gen_samp
We first check the generation with simple sampling of the signals. We see that the shape is not similar.
[6]:
sig3mm.check_gen('lc_gen_samp')
Original vs Synthethic:
mean: 1.8206498796538533 / 1.7996465269173818
std: 1.0533019607954752 / 1.1401532743330116
[6]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
lc_gen_psd_*
[7]:
sig3mm.check_gen('lc_gen_psd_fft')
Original vs Synthethic:
mean: 1.8206498796538533 / 1.8206498796538533
std: 1.0533019607954752 / 1.0533019607954752
/home/docs/checkouts/readthedocs.org/user_builds/mutis/envs/latest/lib/python3.7/site-packages/scipy/signal/spectral.py:1966: UserWarning: nperseg = 256 is greater than input length = 105, using nperseg = 105
.format(nperseg, input_length))
[7]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
We see that the signals have similar shape, and almost identical variability. The statistical distribution is not exactly the same.
We now check using the non-uniform fourier transform. The generation with the FFT was not technically correct since the signal was not evenly sampled in time.
[8]:
sig3mm.check_gen('lc_gen_psd_nft')
Original vs Synthethic:
mean: 1.8206498796538533 / 1.8206498796538535
std: 1.0533019607954752 / 1.053301960795475
[8]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
Still another method can be used to reconstruct the signal; the Lomb-Scargle method to compute the PSD.
[9]:
sig3mm.check_gen('lc_gen_psd_lombscargle')
Original vs Synthethic:
mean: 1.8206498796538533 / 1.8206498796538533
std: 1.0533019607954752 / 1.0530244919105098
[9]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
From these results, we see that the shape for 3mm is better reproduced with ``lc_gen_psd_fft``.
lc_gen_OU
To use the stochastic OU method, first we need to find suitable parameters:
[10]:
sig3mm.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)
/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)
[10]:
{'sigma_est': 2.91683807274248,
'curve_fit': (array([3.02346118, 1.86171047]),
array([0.40840367, 0.11016486])),
'MLE_fit': (2.3531636913848892, 1.9175162668279229),
'th_est1': 10.010292857664975,
'th_est2': 12.861719731625044}
[11]:
sig3mm.check_gen('lc_gen_ou', fgen_params={'mu':1.8, 'sigma':1.01, 'theta':1.96})
Original vs Synthethic:
mean: 1.8206498796538533 / 2.150739942003629
std: 1.0533019607954752 / 1.0867404955082687
[11]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
[12]:
sig3mm.OU_mu, sig3mm.OU_sigma, sig3mm.OU_theta = 1.8, 1.01, 1.96
We see that these value produce synthethic light curves that have similar shape, distribution and variability.
[ ]:
Gamma
[13]:
sigGamma = mutis.Signal(data['gamma']['jyear'][np.isfinite(data['gamma']['CFlux'])], 1e6*data['gamma']['CFlux'][np.isfinite(data['gamma']['CFlux'])], 1e6*data['gamma']['CFluxErr'][np.isfinite(data['gamma']['CFlux'])])
sigGamma.plot()
[13]:
[<matplotlib.lines.Line2D at 0x7ff873948f10>]
lc_gen_psd_*
[14]:
sigGamma.check_gen('lc_gen_psd_fft')
Original vs Synthethic:
mean: 0.11241061716509433 / 0.11241061716509433
std: 0.14036970709552712 / 0.14036970709552712
[14]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
Here we see that lc_gen_psd_fft does not generate signals with similar shapes, specially during flares.
[15]:
sigGamma.check_gen('lc_gen_psd_nft')
Original vs Synthethic:
mean: 0.11241061716509433 / 0.11241061716509433
std: 0.14036970709552712 / 0.14036970709552715
[15]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
lc_gen_OU
[16]:
sigGamma.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)
/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)
[16]:
{'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}
[17]:
sigGamma.check_gen('lc_gen_ou',
fgen_params={'mu':0.21, 'sigma':4.9, 'theta':7.6,
})#'scale':np.std(sigGamma.values), 'loc':np.mean(sigGamma.values)})
Original vs Synthethic:
mean: 0.11241061716509433 / 0.1546240191295745
std: 0.14036970709552712 / 0.2002713001286012
[17]:
{'ax1': <AxesSubplot:title={'center':'light curves'}>,
'ax2': <AxesSubplot:title={'center':'pdf'}>,
'ax3': <AxesSubplot:title={'center':'PSD'}>}
[ ]:
Set and generate synthethics
Now that we have checked which generation methods are the best, we set them and generate the synthethic light curves.
[18]:
%%time
sig3mm.fgen = 'lc_gen_psd_fft'
sig3mm.gen_synth(400)
CPU times: user 121 ms, sys: 4.03 ms, total: 125 ms
Wall time: 125 ms
/home/docs/checkouts/readthedocs.org/user_builds/mutis/envs/latest/lib/python3.7/site-packages/scipy/signal/spectral.py:1966: UserWarning: nperseg = 256 is greater than input length = 105, using nperseg = 105
.format(nperseg, input_length))
[19]:
%%time
sigGamma.fgen = 'lc_gen_ou'
sigGamma.OU_mu, sigGamma.OU_sigma, sigGamma.OU_theta = 0.2, 4.9, 7
sigGamma.gen_synth(400)
CPU times: user 1.42 s, sys: 3.89 ms, total: 1.42 s
Wall time: 1.42 s
Correlation
[20]:
corr3mmGamma = mutis.Correlation(sig3mm, sigGamma, 'welsh')
corr3mmGamma.plot_signals()
[21]:
corr3mmGamma.gen_times(ftimes='uniform', tmin=-200/365, tmax=+200/365, n=50, nbinsmin=11)
corr3mmGamma.plot_times()
plt.xlim([-200/365,+200/365])
[21]:
(-0.547945205479452, 0.547945205479452)
[22]:
corr3mmGamma.samples = 400
[23]:
%%time
corr3mmGamma.gen_corr(uncert=False)
CPU times: user 2.43 s, sys: 0 ns, total: 2.43 s
Wall time: 2.43 s
[24]:
corr3mmGamma.plot_corr(uncert=False)
plt.xlim([-200/365,+200/365])
[24]:
(-0.547945205479452, 0.547945205479452)
[ ]: