Commit 8bc48ed6 authored by Fredrik Bagge Carlson's avatar Fredrik Bagge Carlson
Browse files

Changes to PMMH

parent 826bede6
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)
......
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
......@@ -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.
......
......@@ -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
......
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