In [1]:
push!(LOAD_PATH, "/Users/farr/Documents/Research/GaussianNoiseModeling/PopeFilters/code")
Out[1]:
In [2]:
using Ensemble
using Kalman
using PyCall
using PyPlot
In [3]:
@pyimport seaborn
seaborn.set_style("ticks")
In [4]:
name0 = "EPIC 211000411"
name9 = "EPIC 211098454"
Out[4]:
In [5]:
lc0 = readdlm("epic211000411.txt")
errorbar(lc0[:,1], lc0[:,2], lc0[:,3], fmt=".")
xlabel(L"$t$ ($\mathrm{d}$)")
ylabel(L"$F$ ($\mathrm{cts} \, \mathrm{s}^{-1}$)")
Out[5]:
In [6]:
lc9 = readdlm("epic211098454.txt")
errorbar(lc9[:,1], lc9[:,2], lc9[:,3], fmt=".")
xlabel(L"$t$ ($\mathrm{d}$)")
ylabel(L"$F$ ($\mathrm{cts} \, \mathrm{s}^{-1}$)")
Out[6]:
In [7]:
ns0 = open(deserialize, "epic211000411-CARMA/state-5-4.dat", "r")
ps0, ll0 = EnsembleNest.postsample(ns0)
nps0 = size(ps0, 2)
post0 = Kalman.CARMAKalmanPosterior(lc0[:,1], lc0[:,2], lc0[:,3], 5, 4)
size(ps0)
Out[7]:
In [8]:
plot(ns0.deadlogls[10000:end])
xlabel(L"$N-10000$")
ylabel(L"$\log L$")
Out[8]:
In [8]:
ns9 = open(deserialize, "epic211098454-CARMA//state-5-4.dat", "r")
ps9, ll9 = EnsembleNest.postsample(ns9)
nps9 = size(ps9, 2)
post9 = Kalman.CARMAKalmanPosterior(lc9[:,1], lc9[:,2], lc9[:,3], 5, 4)
size(ps9)
Out[8]:
In [10]:
plot(ns9.deadlogls[10000:end])
xlabel(L"$N-10000$")
ylabel(L"$\log L$")
Out[10]:
In [11]:
fs0 = Kalman.psdfreq(post0)
fs9 = Kalman.psdfreq(post9)
psds0 = Kalman.psd(post0, ps0[:, rand(1:nps0, 1000)], fs0)
psds9 = Kalman.psd(post9, ps9[:, rand(1:nps9, 1000)], fs9)
Out[11]:
Plot the PSDs:
In [12]:
function plot_psd(fs, psd)
nf = size(fs, 1)
ns = size(psd, 2)
m = zeros(nf)
h = zeros(nf)
l = zeros(nf)
for i in 1:nf
m[i] = median(vec(psd[i,:]))
h[i] = quantile(vec(psd[i,:]), 0.84)
l[i] = quantile(vec(psd[i,:]), 0.16)
end
plot(fs, m)
fill_between(fs, h, l, alpha=0.5)
xscale("log")
yscale("log")
xlabel(L"$f$ ($\mathrm{d}^{-1}$)")
ylabel(L"$P(f)$ ($\mathrm{cts}^2 \, \mathrm{s}^{-2} \, \mathrm{d}$)")
end
Out[12]:
The PSD for EPIC 211000411, which has the evolving signal and correspondingly wide modes.
In [13]:
plot_psd(fs0, psds0)
Out[13]:
The PSD for EPIC 211098454, with very sharp modes.
In [14]:
plot_psd(fs9, psds9)
Out[14]:
We want to plot the longest period in the data (lowest frequency) for the two stars, the amplitude ratio of the two peaks that appear, and the period or frequency ratio. For some reason, the seaborn KDE is giving us trouble here, so all the plots are histogram-only.
In [15]:
freqs0 = Kalman.frequencies(post0, ps0)
freqs9 = Kalman.frequencies(post9, ps9)
fmin0 = zeros(size(freqs0,2))
fmax0 = zeros(size(freqs0,2))
fmin9 = zeros(size(freqs9,2))
fmax9 = zeros(size(freqs9,2))
for i in 1:size(fmin0,1)
fs = freqs0[:,i]
fs = fs[fs .> 0]
fmin0[i] = minimum(fs)
fmax0[i] = maximum(fs)
end
for i in 1:size(fmin9,1)
fs = freqs9[:,i]
fs = fs[fs .> 0]
fmin9[i] = minimum(fs)
fmax9[i] = maximum(fs)
end
In [16]:
Pmax0 = 1.0./fmin0
Pmin0 = 1.0./fmax0
Pmax9 = 1.0./fmin9
Pmin9 = 1.0./fmax9
Out[16]:
In [19]:
seaborn.distplot(Pmax0)
xlabel(L"$P_\mathrm{max}$ ($\mathrm{d}$)")
ylabel(L"$p\left( P_\mathrm{max}\right)$")
println("Period of $name0 = $(mean(Pmax0)) +/- $(std(Pmax0)) d")
In [20]:
seaborn.distplot(Pmax9)
xlabel(L"$P_\mathrm{max}$ ($\mathrm{d}$)")
ylabel(L"$p\left( P_\mathrm{max}\right)$")
println("Period of $name9 = $(mean(Pmax9)) +/- $(std(Pmax9)) d")
Note that the "9" period disagrees with the value given in README: it is a factor of two larger. I'm not sure why the smaller period is quoted in the README, but maybe Suzanne reported the stronger of the two modes, no matter which one has larger periods?
How about their ratio?
In [21]:
seaborn.distplot(Pmax0./Pmin0)
xlabel(L"$P_\mathrm{max}/P_\mathrm{min}$")
ylabel(L"$p\left(P_\mathrm{max}/P_\mathrm{min}\right)$")
Out[21]:
In [22]:
seaborn.distplot(Pmax9./Pmin9)
xlabel(L"$P_\mathrm{max}/P_\mathrm{min}$")
ylabel(L"$p\left(P_\mathrm{max}/P_\mathrm{min}\right)$")
Out[22]:
In both systems we seem to have found the first and second harmonics of a signal....
In [23]:
drates0 = Kalman.drates(post0, ps0)
drates9 = Kalman.drates(post9, ps9)
Out[23]:
In [24]:
drate_fund0 = zeros(size(drates0,2))
drate_fund9 = zeros(size(drates9,2))
drate_harm0 = zeros(size(drates0,2))
drate_harm9 = zeros(size(drates9,2))
for i in 1:size(drates0, 2)
inds = sortperm(freqs0[:,i])
drate_fund0[i] = drates0[inds[end-1],i] # Smallest (positive) freq decay rate
drate_harm0[i] = drates0[inds[end],i] # Largest freq decay rate
end
for i in 1:size(drates9, 2)
inds = sortperm(freqs9[:,i])
drate_fund9[i] = drates9[inds[end-1],i]
drate_harm9[i] = drates9[inds[end],i]
end
Timescales assocaited with the fundamental and harmonic mode for EPIC 211000411. The prior on the timescale extends only out to $T_\mathrm{obs}$, so the really sharp (i.e. long-timescale) modes have a cutoff due to the prior.
You can see that the timescale for the fundamental mode is permitted to be quite short---i.e. that the mode is evolving, as is evident in the lightcurve.
In [25]:
seaborn.distplot(-1.0./drate_fund0, label="Fundamental")
seaborn.distplot(-1.0./drate_harm0, label="Harmonic")
xlabel(L"$\tau$ ($\mathrm{d}$)")
ylabel(L"$p\left(\tau\right)$")
legend(loc="upper left")
Out[25]:
Timescales associated with the fundamential and harmonic mode of EPIC 211098454
In [27]:
seaborn.distplot(-1.0./drate_fund9, label="Fundamental")
seaborn.distplot(-1.0./drate_harm9, label="Harmonic")
xlabel(L"$\tau$ ($\mathrm{d}$)")
ylabel(L"$p\left(\tau\right)$")
legend(loc="upper left")
Out[27]:
TODO: Amplitude ratios!
These models allow for three modes.
In [9]:
ns976 = open(deserialize, "epic211098454-CARMA/state-7-6.dat", "r")
ns076 = open(deserialize, "epic211000411-CARMA/state-7-6.dat", "r")
ps076, ll076 = EnsembleNest.postsample(ns076)
ps976, ll976 = EnsembleNest.postsample(ns976)
post076 = Kalman.CARMAKalmanPosterior(lc0[:,1], lc0[:,2], lc0[:,3], 7, 6)
post976 = Kalman.CARMAKalmanPosterior(lc9[:,1], lc9[:,2], lc9[:,3], 7, 6)
nps076 = size(ps076, 2)
nps976 = size(ps976, 2)
Out[9]:
Are the fits converged?
In [17]:
plot(ns976.deadlogls[10000:end])
Out[17]:
In [18]:
plot(ns076.deadlogls[10000:end])
Out[18]:
In [29]:
psds076 = Kalman.psd(post076, ps076[:, rand(1:nps076, 1000)], fs0)
psds976 = Kalman.psd(post976, ps976[:, rand(1:nps976, 1000)], fs9)
Out[29]:
Below in the $(7,6)$ estimate of the PSD; you can see that there is no evidence for another mode in the data, though there seems to be more low-frequency power than appeared in the $(5,4)$ fit. There is probably low-frequency fluctuations (maybe associated to imperfect systematic removal?) present in the data at lower power than the modes; the $(5,4)$ model is essentially tapped out after accounting for the two modes (4 out of 5 AR parameters) and the falloff at high frequency (the remaining AR parameter).
In [30]:
plot_psd(fs0, psds076)
plot_psd(fs9, psds976)
legend([name0, name9], loc="upper right")
Out[30]:
What are the evidences for $(7,6)$ versus $(5,4)$?
In [10]:
println("$(name0): ln(Z_76/Z_54) = $(EnsembleNest.logZ(ns076)-EnsembleNest.logZ(ns0))")
println("$(name9): ln(Z_76/Z_54) = $(EnsembleNest.logZ(ns976)-EnsembleNest.logZ(ns9))")
How about the DIC:
In [12]:
println("$(name0): DIC(76) - DIC(54) = $(EnsembleNest.dic(ns076) - EnsembleNest.dic(ns0))")
println("$(name9): DIC(76) - DIC(54) = $(EnsembleNest.dic(ns976)-EnsembleNest.dic(ns9))")
So, the DIC is going in the opposite direction to the evidence in EPIC 211000411 (smaller DIC is better than larger), while both measures prefer the $(5,4)$ model for EPIC 211098454.
In [ ]: