Commit e6c4b901 authored by Fredrik Bagge Carlson's avatar Fredrik Bagge Carlson
Browse files

Updates to particle filter implementations

parent 9c47a071
...@@ -28,7 +28,7 @@ f(xn,x,t::Int64) = begin ...@@ -28,7 +28,7 @@ f(xn,x,t::Int64) = begin
end end
end end
function f_sample(x::Vector, t::Int64) function f_sample(x, t::Int64)
c = 8*cos(1.2*(t-1)) c = 8*cos(1.2*(t-1))
@inbounds for i = 1:length(x) @inbounds for i = 1:length(x)
x[i] = 0.5*x[i] + 25*x[i]./(1+x[i]^2) + c + σw*randn() x[i] = 0.5*x[i] + 25*x[i]./(1+x[i]^2) + c + σw*randn()
...@@ -45,19 +45,19 @@ function f_density(xb, xp, t) ...@@ -45,19 +45,19 @@ function f_density(xb, xp, t)
end end
const den = 0.5/σv^2 const den = -0.5/σv^2
function g_density(y::Float64, x::Vector{Float64}, w::Vector{Float64}) function g_density(y::Float64, x, w)
@inbounds for i = 1:length(w) @inbounds for i = 1:length(w)
w[i] -= (den * (y-0.05x[i]^2)^2 - s2piσv) w[i] += (den * (y-0.05x[i]^2)^2 - s2piσv)
end end
end end
## Parameter estimation ## Parameter estimation
function f_theta_sample(x::Vector, t, theta) function f_theta_sample(x,x1, t, theta)
c = theta[3]*cos(1.2*(t-1)) c = theta[3]*cos(1.2*(t-1))
@inbounds for i = 1:length(x) @inbounds for i = 1:length(x)
x[i] = theta[1]*x[i] + theta[2]*x[i]./(1+x[i]^2) + c + σw*randn() x[i] = theta[1]*x1[i] + theta[2]*x1[i]./(1+x1[i]^2) + c + σw*randn()
end end
return x return x
end end
...@@ -72,10 +72,10 @@ function f_density_theta(xb, xp, t, theta) ...@@ -72,10 +72,10 @@ function f_density_theta(xb, xp, t, theta)
end end
function g_density_theta(y::Float64,x::Vector{Float64},w, theta) function g_density_theta(y::Float64,x, w, theta)
# w = Array(Float64,size(x)) # w = Array(Float64,size(x))
@inbounds for i = 1:length(w) @inbounds for i = 1:length(w)
w[i] = -(den * (y-0.05x[i]^2)^2) w[i] = den * (y-0.05x[i]^2)^2 - s2piσv
end end
return w return w
end end
...@@ -133,14 +133,14 @@ function test_pf() ...@@ -133,14 +133,14 @@ function test_pf()
return RMSE, RMSE_FFBSi return RMSE, RMSE_FFBSi
end end
using Winston
function test_PMMH() function test_PMMH()
particle_count = [100] particle_count = [100]
time_steps = [200] time_steps = [200]
R = 5_000 R = 30_000
for (Ti,T) in enumerate(time_steps) for (Ti,T) in enumerate(time_steps)
for (Ni, N) in enumerate(particle_count) for (Ni, N) in enumerate(particle_count)
montecarlo_runs = 4# maximum(particle_count)*maximum(time_steps) / T / N*10 montecarlo_runs = 8# maximum(particle_count)*maximum(time_steps) / T / N*10
x = Array(Float64,T) x = Array(Float64,T)
y = Array(Float64,T) y = Array(Float64,T)
theta_PMMH = SharedArray(Float64,(size(theta0,1),R,montecarlo_runs)) theta_PMMH = SharedArray(Float64,(size(theta0,1),R,montecarlo_runs))
...@@ -152,13 +152,18 @@ function test_PMMH() ...@@ -152,13 +152,18 @@ function test_PMMH()
y[t+1] = 0.05x[t+1]^2 + σv*randn() y[t+1] = 0.05x[t+1]^2 + σv*randn()
end # t end # t
theta0i = draw_theta0() theta0i = draw_theta0()
theta_PMMH[:,:,mc_iter] = pf_PMMH(y, N, R, pf_bootstrap_nn, f_theta_sample, draw_theta_random, f_density_theta, g_density_theta, σw0, theta0i) theta_PMMH[:,:,mc_iter] = pf_PMMH(y, N, R, pf_bootstrap_nn, f_theta_sample, draw_theta_random, f_density_theta, g_density_theta, theta0i)
end # MC end # MC
newplot(theta_PMMH[:,:,1]'); hold(true) newplot(theta_PMMH[:,:,1]'); hold(true)
for mc_iter = 2:montecarlo_runs for mc_iter = 2:montecarlo_runs
plot(theta_PMMH[:,:,mc_iter]') plot(theta_PMMH[:,:,mc_iter]')
end end
hold(false); title("PMMH") hold(false); title("PMMH")
th = Array(theta_PMMH)
# pp = Winston.figure(); Winston.plothist(th[1,:,:][:],500); Winston.hold(true)
# Winston.plothist(th[2,:,:][:],500)
# Winston.plothist(th[3,:,:][:],500)
# Winston.hold(false); Winston.display(pp)
@show N @show N
end # N end # N
@show T @show T
......
...@@ -12,7 +12,7 @@ function pf_PMMH(y, N, R, pf, f_theta, draw_theta, f_density_theta, g_density_th ...@@ -12,7 +12,7 @@ function pf_PMMH(y, N, R, pf, f_theta, draw_theta, f_density_theta, g_density_th
theta[:,1] = theta0 theta[:,1] = theta0
g_density(y_t, x_t,w) = g_density_theta(y_t, x_t,w, theta0) g_density(y_t, x_t,w) = g_density_theta(y_t, x_t,w, theta0)
f(x,t) = f_theta(x,t,theta0) f(x,x1,t) = f_theta(x,x1,t,theta0)
pf(xp, w, wnn, y, N, g_density, f) pf(xp, w, wnn, y, N, g_density, f)
Z[1] = Zf(wnn) Z[1] = Zf(wnn)
...@@ -24,7 +24,7 @@ function pf_PMMH(y, N, R, pf, f_theta, draw_theta, f_density_theta, g_density_th ...@@ -24,7 +24,7 @@ function pf_PMMH(y, N, R, pf, f_theta, draw_theta, f_density_theta, g_density_th
breakTime = floor(Int64, R / breakPoints) breakTime = floor(Int64, R / breakPoints)
particleIncrease = 10*N particleIncrease = 10*N
thetaP = draw_theta(theta0) thetaP = draw_theta(theta0)
f(x,t) = f_theta(x,t,thetaP) f(x,x1,t) = f_theta(x,x1,t,thetaP)
g_density(y_t, x_t,w) = g_density_theta(y_t, x_t,w, thetaP) g_density(y_t, x_t,w) = g_density_theta(y_t, x_t,w, thetaP)
for r = 2:R for r = 2:R
......
...@@ -7,7 +7,7 @@ function pf_bootstrap!(xp,w, y, N, g_density, f) ...@@ -7,7 +7,7 @@ function pf_bootstrap!(xp,w, y, N, g_density, f)
wT = slice(w,:,1) wT = slice(w,:,1)
xT = slice(xp,:,1) xT = slice(xp,:,1)
fill!(wT,-log(N)) fill!(wT,-log(N))
g(y[1],xT, wT) g_density(y[1],xT, wT)
wT -= logsumexp(wT) wT -= logsumexp(wT)
j = Array(Int64,N) j = Array(Int64,N)
# tw = Float64(N) # tw = Float64(N)
...@@ -27,7 +27,7 @@ function pf_bootstrap!(xp,w, y, N, g_density, f) ...@@ -27,7 +27,7 @@ function pf_bootstrap!(xp,w, y, N, g_density, f)
end end
# Time update # Time update
g(y[t],xT, wT) g_density(y[t],xT, wT)
# Normalize weights # Normalize weights
wT -= logsumexp(wT) wT -= logsumexp(wT)
...@@ -45,7 +45,7 @@ function pf_bootstrap(y, N, g_density, f) ...@@ -45,7 +45,7 @@ function pf_bootstrap(y, N, g_density, f)
wT = slice(w,:,1) wT = slice(w,:,1)
xT = slice(xp,:,1) xT = slice(xp,:,1)
fill!(wT,-log(N)) fill!(wT,-log(N))
g(y[1],xT, wT) g_density(y[1],xT, wT)
wT -= logsumexp(wT) wT -= logsumexp(wT)
j = Array(Int64,N) j = Array(Int64,N)
# tw = Float64(N) # tw = Float64(N)
...@@ -65,7 +65,7 @@ function pf_bootstrap(y, N, g_density, f) ...@@ -65,7 +65,7 @@ function pf_bootstrap(y, N, g_density, f)
end end
# Time update # Time update
g(y[t],xT, wT) g_density(y[t],xT, wT)
# Normalize weights # Normalize weights
wT -= logsumexp(wT) wT -= logsumexp(wT)
...@@ -84,8 +84,8 @@ function pf_bootstrap_nn(xp, w, wnn, y, N, g_density, f) ...@@ -84,8 +84,8 @@ function pf_bootstrap_nn(xp, w, wnn, y, N, g_density, f)
wT = slice(w,:,1) wT = slice(w,:,1)
xT = slice(xp,:,1) xT = slice(xp,:,1)
fill!(wT,-log(N)) fill!(wT,-log(N))
g(y[1],xT, wT) g_density(y[1],xT, wT)
wnn[:,t] = wT wnn[:,1] = wT
wT -= logsumexp(wT) wT -= logsumexp(wT)
j = Array(Int64,N) j = Array(Int64,N)
# tw = Float64(N) # tw = Float64(N)
...@@ -178,7 +178,7 @@ function pf_CSMC!(xp,w, y, N, g_density, f, xc ) ...@@ -178,7 +178,7 @@ function pf_CSMC!(xp,w, y, N, g_density, f, xc )
wT = slice(w,:,1) wT = slice(w,:,1)
xT = slice(xp,:,1) xT = slice(xp,:,1)
fill!(wT,-log(N)) fill!(wT,-log(N))
g(y[1],xT, wT) g_density(y[1],xT, wT)
wT -= logsumexp(wT) wT -= logsumexp(wT)
j = Array(Int64,N) j = Array(Int64,N)
# tw = Float64(N) # tw = Float64(N)
...@@ -198,7 +198,7 @@ function pf_CSMC!(xp,w, y, N, g_density, f, xc ) ...@@ -198,7 +198,7 @@ function pf_CSMC!(xp,w, y, N, g_density, f, xc )
end end
xp[N,t] = xc[t] xp[N,t] = xc[t]
# Time update # Time update
g(y[t],xT, wT) g_density(y[t],xT, wT)
# Normalize weights # Normalize weights
wT -= logsumexp(wT) wT -= logsumexp(wT)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment