In [1]:
push!(LOAD_PATH, ".")
Out[1]:
In [2]:
using DataFrames
using Gadfly
using Kalman
using Base.Test:@test_approx_eq_eps
Generate some non-uniformly spaced times, and some observational errors:
In [3]:
ts = collect(1:3000) + 0.5*rand(3000) - 0.25
dys = 0.1 + 0.1*rand(3000)
nothing
In [4]:
dt = minimum(diff(ts))
T = (ts[end]-ts[1])
nothing
In [5]:
function randtscale(dt, T)
exp(log(dt) + rand()*log(T/dt))
end
Out[5]:
In [6]:
mu = cbrt(pi)
sigma = float(pi)
x = -randtscale(dt, T)
y = randtscale(dt, T)
arroots = [-1.0/randtscale(dt,T), 1.0/(x+1im*y), 1.0/(x-1im*y)]
x = randtscale(dt,T)
y = randtscale(dt,T)
maroots = [1.0/(-x+1im*y), 1.0/(-x-1im*y)]
nothing
In [7]:
filt = Kalman.CARMAKalmanFilter(mu, sigma, arroots, maroots)
sfilt = Kalman.CARMASqrtKalmanFilter(mu, sigma, arroots, maroots)
nothing
In [8]:
ys = Kalman.carmagenerate(ts, dys, mu, sigma, arroots, maroots)
plot(x=ts, y=ys, Geom.line)
Out[8]:
Let's look at the equivalent generated from the filter (these won't match exactly, of course, but they will look similar, hopefully):
In [15]:
ysg = Kalman.generate(filt, ts, dys)
println("Mean of ysg = $(mean(ysg)), mu = $mu")
println("Std of ysg = $(std(ysg)), sigma = $sigma")
plot(x=ts, y=ysg, Geom.line)
Out[15]:
In [9]:
Kalman.raw_carma_log_likelihood(ts, ys, dys, mu, sigma, arroots, maroots)
Out[9]:
In [10]:
Kalman.log_likelihood(sfilt, ts, ys, dys)
Out[10]:
In [11]:
Kalman.log_likelihood(filt, ts, ys, dys)
Out[11]:
In [12]:
@time for i in 1:1000; Kalman.log_likelihood(filt, ts, ys, dys) end
In [13]:
Profile.clear()
@profile for i in 1:1000; Kalman.log_likelihood(filt, ts, ys, dys) end
Profile.print()
In [ ]: