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 corner
@pyimport seaborn as sns
sns.set_style("ticks")
In [4]:
data = readdlm("../epic210969800.txt")
errorbar(data[:,1], data[:,2], data[:,3], fmt=".")
xlabel(L"$t$ ($\mathrm{d}$)")
ylabel(L"$F$ ($\mathrm{cts} \, \mathrm{s}^{-1}$)")
Out[4]:
In [5]:
post = Kalman.CARMAKalmanPosterior(data[:,1], data[:,2], data[:,3], 3, 2)
function logpost(x)
lp = Kalman.log_prior(post, x)
if lp == -Inf
lp
else
lp + Kalman.log_likelihood(post, x)
end
end
Out[5]:
In [6]:
pts = Kalman.init(post, 128)
lnprobs = EnsembleSampler.lnprobs(pts, logpost)
Out[6]:
In [7]:
ps_save = nothing
lnps_save = nothing
In [8]:
function neff_cb(ps, lnps, n, thin)
println("Iterated by $(n)")
println(" (max, mean, var)(log(L)) = $((maximum(lnps), mean(lnps), var(lnps)))")
println(" max ACL = $(maximum(Acor.acl(ps)))")
ps_save = ps
lnps_save = lnps
end
@time ens, enslnps = EnsembleSampler.run_to_neff(pts, lnprobs, logpost, 16, callback=neff_cb)
Out[8]:
The above run involved $\mathcal{O}(3000)$ ensemble steps, which is $\mathcal{O}(300000)$ likelihood evaluations. The equivalent nested sampling run retired $33792$ points, with each one taking $\mathcal{O}(100)$ likelihood evaluations, which is a factor of ten or more longer!
In [9]:
Acor.gelman_rubin_rs(ens)
Out[9]:
In [10]:
for i in 1:128
plot(vec(enslnps[i,:]))
end
In [14]:
for i in 1:8
subplot(3,3,i)
plot(vec(mean(ens[i,:,:], 2)))
for j in 1:10
k = rand(1:128)
plot(vec(ens[i,k,:]), alpha=0.1)
end
end
In [19]:
r, sr = Kalman.residuals(post, ens[:,1,128])
plot(data[:,1], r./sr)
figure()
sns.kdeplot(r./sr)
xs = collect(-3:0.01:3)
plot(xs, 1.0/sqrt(2.0*pi)*exp(-0.5*xs.*xs))
Out[19]:
In [21]:
ts = sort(vcat(data[:,1], collect(linspace(minimum(data[:,1]), maximum(data[:,1]), 1000))))
ys, vys = Kalman.predict(post, ens[:,1,128], ts)
plot(ts, ys)
plot(data[:,1], data[:,2], ".k")
fill_between(ts, ys+sqrt(vys), ys-sqrt(vys), alpha=0.5)
Out[21]:
In [25]:
fs = Kalman.psdfreq(post)
psds = Kalman.psd(post, reshape(ens, (size(ens,1), size(ens,2)*size(ens,3))), fs)
mid = zeros(size(psds,1))
low = zeros(size(psds,1))
high = zeros(size(psds,1))
for i in 1:size(psds,1)
mid[i] = median(vec(psds[i,:]))
low[i] = quantile(vec(psds[i,:]), 0.16)
high[i] = quantile(vec(psds[i,:]), 0.84)
end
plot(fs, mid)
fill_between(fs, high, low, alpha=0.5)
xscale("log")
yscale("log")
In [ ]: