From 8bc48ed63113601144c3351909e94c809ebbc6e4 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <fredrikb@control.lth.se>
Date: Sun, 13 Sep 2015 23:29:23 +0200
Subject: [PATCH] Changes to PMMH
---
.../particle_smoother_RS-FFBSi.jl | 1 +
src/particle_filters/pf_PMMH.jl | 68 +++++++++++++++++++
src/particle_filters/resample.jl | 20 +++++-
src/sobol.jl | 6 +-
4 files changed, 89 insertions(+), 6 deletions(-)
create mode 100644 src/particle_filters/pf_PMMH.jl
diff --git a/src/particle_filters/particle_smoother_RS-FFBSi.jl b/src/particle_filters/particle_smoother_RS-FFBSi.jl
index 2f725a8..4e5a464 100644
--- a/src/particle_filters/particle_smoother_RS-FFBSi.jl
+++ b/src/particle_filters/particle_smoother_RS-FFBSi.jl
@@ -1,4 +1,5 @@
function pf_smoother_FFBSi( w,xp, M, f_density )
+ error("Not yet implemented")
(N,T) = size(w)
xb = Array(Float64,M,T)
wexp = exp(w)
diff --git a/src/particle_filters/pf_PMMH.jl b/src/particle_filters/pf_PMMH.jl
new file mode 100644
index 0000000..b87b5f0
--- /dev/null
+++ b/src/particle_filters/pf_PMMH.jl
@@ -0,0 +1,68 @@
+using Debug
+function pf_PMMH(y, N, R, pf, f_theta, draw_theta, f_density_theta, g_density_theta, sigmaW0, theta0 )
+
+ n_params = length(theta0)
+ T = length(y)
+ xp = Array(Float64,(N,T))
+ w = Array(Float64,(N,T))
+ wnn = Array(Float64,(N,T))
+ Zf(wnn) = sum(log(mean(exp(wnn),1)))
+ Z = Array(Float64,R)
+
+ theta = Array(Float64,n_params,R)
+ 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)
+ pf(xp, w, wnn, y, N, g_density, f, sigmaW0)
+
+ Z[1] = Zf(wnn)
+
+ breakPoints = 10
+
+ # noiseMultiplier = 0.1^(1/R)
+ Nstart = N
+ breakTime = floor(Int64, R / breakPoints)
+ particleIncrease = 10*N
+ thetaP = draw_theta(theta0)
+ f(x,t) = f_theta(x,t,thetaP)
+ g_density(y_t, x_t,w) = g_density_theta(y_t, x_t,w, thetaP)
+
+ for r = 2:R
+ thetaP = draw_theta(theta[:,r-1])
+
+ # f_density(x_t1, x_t, t) = f_density_theta(x_t1, x_t, t, thetaP)
+
+
+ pf(xp, w, wnn, y, N, g_density, f, sigmaW0)
+ Zp = Zf(wnn)
+# @show exp(Zp-Z[r-1])
+ Iaccept = rand() < exp(Zp-Z[r-1])
+
+ if Iaccept
+ theta[:,r] = thetaP
+ Z[r] = Zp
+ else
+ theta[:,r] = theta[:,r-1]
+ Z[r] = Z[r-1]
+ end
+
+ if (r % breakTime == 0)
+ display(r/R)
+# N = round(Int64, particleIncrease*(r/R)^3 + Nstart )
+# xp = Array(Float64,(N,T))
+# w = Array(Float64,(N,T))
+# wnn = Array(Float64,(N,T))
+ end
+
+
+ end
+
+ return theta
+end
+
+
+
+
+
+
diff --git a/src/particle_filters/resample.jl b/src/particle_filters/resample.jl
index a72fa33..293b90c 100644
--- a/src/particle_filters/resample.jl
+++ b/src/particle_filters/resample.jl
@@ -22,7 +22,7 @@ function resample_systematic!(j,w,M)
s = collect((rand()/M+0):1/M:bins[end])
bo = 1
for i = 1:M
- @inbounds for b = bo:N
+ for b = bo:N
if s[i] < bins[b]
j[i] = b
bo = b
@@ -51,6 +51,24 @@ function resample_systematic_exp(w,M)
return j
end
+function resample_systematic_exp!(j,w)
+ N = size(w,1)
+ M = size(j,1)
+ bins = cumsum(w) # devec
+ s = collect((rand()/M+0):1/M:bins[end])
+ bo = 1
+ for i = 1:M
+ for b = bo:N
+ if s[i] <= bins[b]
+ j[i] = b
+ bo = b
+ break
+ end
+ end
+ end
+ return j
+end
+
# """
# There is probably lots of room for improvement here. All bins need not be formed in the beginning.
# One only has to keep 1 values, the current upper limit, no array needed.
diff --git a/src/sobol.jl b/src/sobol.jl
index 64b455d..5ab9290 100644
--- a/src/sobol.jl
+++ b/src/sobol.jl
@@ -44,9 +44,6 @@ function sobol(X, nextseed, MeM)
return X, nextseed, MeM
end
-## Functions below are inner functions, I have not determined if this is needed, the fastest method or anything else
-# ------------------------------------------------------------------------
-
function getNewSobolVector(dim, seed, MeM)
seed::Int64 = max(floor(seed),0)
@@ -95,7 +92,6 @@ function getNewSobolVector(dim, seed, MeM)
return qrv, nextseed, MeM
end
-# -------------------------------------------------------------------------
function smallest0bit(b)
i = 0
@@ -111,7 +107,7 @@ function smallest0bit(b)
return i
end
-# ----------------------------------------------------------------------
+
function InitSobol(dim::Int64)
MeM = Mem(dim)
MeM.seed_save = -1
--
GitLab