diff --git a/src/particle_filters/particle_filter_test.jl b/src/particle_filters/particle_filter_test.jl
index 3698ee7d3bc268d2f80376e576218626b09e1c42..c20263ed2b574a1996e08fa903d60340d83efe27 100644
--- a/src/particle_filters/particle_filter_test.jl
+++ b/src/particle_filters/particle_filter_test.jl
@@ -28,7 +28,7 @@ f(xn,x,t::Int64) = begin
end
end
-function f_sample(x::Vector, t::Int64)
+function f_sample(x, t::Int64)
c = 8*cos(1.2*(t-1))
@inbounds for i = 1:length(x)
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)
end
-const den = 0.5/σv^2
-function g_density(y::Float64, x::Vector{Float64}, w::Vector{Float64})
+const den = -0.5/σv^2
+function g_density(y::Float64, x, 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
## 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))
@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
return x
end
@@ -72,10 +72,10 @@ function f_density_theta(xb, xp, t, theta)
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))
@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
return w
end
@@ -133,14 +133,14 @@ function test_pf()
return RMSE, RMSE_FFBSi
end
-using Winston
+
function test_PMMH()
particle_count = [100]
time_steps = [200]
- R = 5_000
+ R = 30_000
for (Ti,T) in enumerate(time_steps)
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)
y = Array(Float64,T)
theta_PMMH = SharedArray(Float64,(size(theta0,1),R,montecarlo_runs))
@@ -152,13 +152,18 @@ function test_PMMH()
y[t+1] = 0.05x[t+1]^2 + σv*randn()
end # t
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
newplot(theta_PMMH[:,:,1]'); hold(true)
for mc_iter = 2:montecarlo_runs
plot(theta_PMMH[:,:,mc_iter]')
end
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
end # N
@show T
diff --git a/src/particle_filters/pf_PMMH.jl b/src/particle_filters/pf_PMMH.jl
index 6f3dd0e636585a90f8e7aac19a74526b3b880491..f97e3ff96f2b3d86543ce82010416ce632e81f79 100644
--- a/src/particle_filters/pf_PMMH.jl
+++ b/src/particle_filters/pf_PMMH.jl
@@ -12,7 +12,7 @@ function pf_PMMH(y, N, R, pf, f_theta, draw_theta, f_density_theta, g_density_th
theta[:,1] = 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)
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
breakTime = floor(Int64, R / breakPoints)
particleIncrease = 10*N
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)
for r = 2:R
diff --git a/src/particle_filters/pf_bootstrap.jl b/src/particle_filters/pf_bootstrap.jl
index 0cb4a6126ee2fff7f3ae82a67e2da5d9fe70b128..931280d76efed8fb248d68a2f5b190016d27171b 100644
--- a/src/particle_filters/pf_bootstrap.jl
+++ b/src/particle_filters/pf_bootstrap.jl
@@ -7,7 +7,7 @@ function pf_bootstrap!(xp,w, y, N, g_density, f)
wT = slice(w,:,1)
xT = slice(xp,:,1)
fill!(wT,-log(N))
- g(y[1],xT, wT)
+ g_density(y[1],xT, wT)
wT -= logsumexp(wT)
j = Array(Int64,N)
# tw = Float64(N)
@@ -27,7 +27,7 @@ function pf_bootstrap!(xp,w, y, N, g_density, f)
end
# Time update
- g(y[t],xT, wT)
+ g_density(y[t],xT, wT)
# Normalize weights
wT -= logsumexp(wT)
@@ -45,7 +45,7 @@ function pf_bootstrap(y, N, g_density, f)
wT = slice(w,:,1)
xT = slice(xp,:,1)
fill!(wT,-log(N))
- g(y[1],xT, wT)
+ g_density(y[1],xT, wT)
wT -= logsumexp(wT)
j = Array(Int64,N)
# tw = Float64(N)
@@ -65,7 +65,7 @@ function pf_bootstrap(y, N, g_density, f)
end
# Time update
- g(y[t],xT, wT)
+ g_density(y[t],xT, wT)
# Normalize weights
wT -= logsumexp(wT)
@@ -84,8 +84,8 @@ function pf_bootstrap_nn(xp, w, wnn, y, N, g_density, f)
wT = slice(w,:,1)
xT = slice(xp,:,1)
fill!(wT,-log(N))
- g(y[1],xT, wT)
- wnn[:,t] = wT
+ g_density(y[1],xT, wT)
+ wnn[:,1] = wT
wT -= logsumexp(wT)
j = Array(Int64,N)
# tw = Float64(N)
@@ -178,7 +178,7 @@ function pf_CSMC!(xp,w, y, N, g_density, f, xc )
wT = slice(w,:,1)
xT = slice(xp,:,1)
fill!(wT,-log(N))
- g(y[1],xT, wT)
+ g_density(y[1],xT, wT)
wT -= logsumexp(wT)
j = Array(Int64,N)
# tw = Float64(N)
@@ -198,7 +198,7 @@ function pf_CSMC!(xp,w, y, N, g_density, f, xc )
end
xp[N,t] = xc[t]
# Time update
- g(y[t],xT, wT)
+ g_density(y[t],xT, wT)
# Normalize weights
wT -= logsumexp(wT)