- First, we load numpy as the clag package. pylab is for plotting
In [10]:
import numpy as np
import clag
%pylab inline
- Now, we read the first light curve. e.g. 1367A.dat, and plot it
In [2]:
dt = 0.01
t1, l1, l1e = np.loadtxt('1367A.dat').T
errorbar(t1, l1, yerr=l1e, fmt='o')
Out[2]:
- Here, we define the frequency bins in fqL. Here, I use the ones that Ed sent.
A general rules for fqL is as follows:
- define $f_1$, $f_2$ as the two extreme frequencies allowed by the data. i.e. $f_{1} = 1/T$ with $T$ being the length of observation in time units, and $f_{2} = 0.5/\Delta t$
- The frequency limits where the psd/lag can be constrained are about ~$0.9f_1 - 0.5f_2$. The 0.9 factor doesn't depend on the data much, but values in the range ~[0.9-1.1] are ok. The factor in front of $f_2$ depends on the quality of the data, for low qualily data, use ~0.1--0.2, and for high quality data, increase it up to $0.9--1$.
- Always include two dummy bins at the low and high frequencies and ignore them. The first and last bins are generally biased, So I suggest using the first bin as $0.5f_1 - 0.9 f_1$ (or whatever value you use instead of $0.9 f_1$, see second point above), and the last bin should be $0.5f_2-2*f_2$ (or whatever value instead of $0.5f_2$, see second point above). So the frequency bins should be something like: $[0.5f_1, 0.9f_1, ..., 0.5f_2, 2f_2]$, the bins in between can be devided as desired.
fqd is the bin center
In [5]:
fqL = np.array([0.0049999999, 0.018619375, 0.044733049, 0.069336227, 0.10747115, 0.16658029,
0.25819945, 0.40020915, 0.62032418])
nfq = len(fqL) - 1
fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.)
- We create a psd class of type psd10r (i.e. psd in log10 units using the rms normalization, see tutorial of clag for more). Other paramters are given as a list of: time, light curve, error. If you want to fit more than one light curve together, give them as elements of that list.
- p1 is an array of the starting parameters
- We then call clag.optimize to find the best fit. It takes as parameters, the class of psd or lag, the starting paramters, plus some other parameters whose default values work in most cases. The function returns a tuple of two arrays: the best fit values, and the estimated errors from the Fisher matrix (i.e. approximate errors).
- As the fit starts, progress is printed, where for every line, we print: num_of_iteration, max change in parameter values, maximum gradient, change in loglikelihood -- loglikelihood -- list of parameters.
- The fit stops when either of the three numbers after iter_number is close to zero.
- At the end, the function prints the best fit values, their estimated error and the gradient of the loglikelihood at the those values. For the fit to be valid, the gradient values need to be all close to zero
In [6]:
P1 = clag.clag('psd10r', [t1], [l1], [l1e], dt, fqL)
p1 = np.ones(nfq)
p1, p1e = clag.optimize(P1, p1)
- We then call clag.errors to obtain the errors on the parameters. This steps over the parameters, allowing others to change, until the loglikelihood function changes by 1. This gives the 1-sigma error.
- The errors can also be obtained using mcmc, though not directly implemeted in clag code, one can use other mcmc packages like emcee, and call the p1.loglikelihood function directly.
In [14]:
p1, p1e = clag.errors(P1, p1, p1e)
- Next we plot the results
In [76]:
xscale('log'); ylim(-4,2)
errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10)
Out[76]:
In [ ]:
- Now, we load the second light curve and repeat the process for the second psd
In [15]:
t2, l2, l2e = np.loadtxt('3465A.dat').T
errorbar(t1, l1, yerr=l1e, fmt='o')
errorbar(t2, l2, yerr=l2e, fmt='o')
Out[15]:
In [16]:
P2 = clag.clag('psd10r', [t2], [l2], [l2e], dt, fqL)
p2 = np.ones(nfq)
p2, p2e = clag.optimize(P2, p2)
In [79]:
xscale('log'); ylim(-6,2)
errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10)
errorbar(fqd, p2, yerr=p2e, fmt='o', ms=10)
Out[79]:
In [17]:
p2, p2e = clag.errors(P2, p2, p2e)
In [ ]:
- Next we creat a lag object similar to the psd case. Now we use both light curves and the psd values obtained earlier.
- We use cxd10r, becasue we used psd10r above, if that changes, this should change too. Though generally, there is no need to change these.
- The light curves, time and error are passed as lists of lists. The outer list contains how many light curve sets we want to fit together, and the inner list contrain the first and second light curves whose lag is to be calculated
- p is an array of starting parameters. A start like that should using works fine. Sometimes this need to be tweaked to make sure the fit converges.
In [52]:
Cx = clag.clag('cxd10r', [[t1,t2]], [[l1,l2]], [[l1e,l2e]], dt, fqL, p1, p2)
p = np.concatenate( ((p1+p2)*0.5-0.3,p1*0+0.1) ) # a good starting point generally
p, pe = clag.optimize(Cx, p)
- We get the errors for the cross spectrum and phase lags
In [96]:
%autoreload
p, pe = clag.errors(Cx, p, pe)
In [99]:
phi, phie = p[nfq:], pe[nfq:]
lag, lage = phi/(2*np.pi*fqd), phie/(2*np.pi*fqd)
cx, cxe = p[:nfq], pe[:nfq]
In [100]:
xscale('log'); ylim(-10,10)
errorbar(fqd, lag, yerr=lage, fmt='o', ms=10)
Out[100]:
In [ ]: