Commit 826bede6 authored by Fredrik Bagge Carlson's avatar Fredrik Bagge Carlson
Browse files

Added some particle filter algorithms, not the optimized filter though

parent 70a37ce3
function pf_smoother_FFBSi( w,xp, M, f_density )
(N,T) = size(w)
xb = Array(Float64,M,T)
wexp = exp(w)
j = resample_systematic_exp(wexp[:,T],M)
xb[:,T] = xp[j, T]
wb = Array(Float64,N)
for t = T-1:-1:1
for m = 1:M
for n = 1:N
wb[n] = wexp[n,t]*f_density(xb[m,t+1],xp[n,t],t)+eps()
end
i = draw_one_categorical(wb)
xb[m,t] = xp[i, t]
end
end
return xb
end
function pf_smoother_FFBSi( w,xp, M, f_density )
(N,T) = size(w)
xb = Array(Float64,M,T)
wexp = exp(w)
j = resample_systematic_exp(wexp[:,T],M)
xb[:,T] = xp[j, T]
wb = Array(Float64,N)
for t = T-1:-1:1
for m = 1:M
for n = 1:N
wb[n] = wexp[n,t]*f_density(xb[m,t+1],xp[n,t],t)+eps()
end
i = draw_one_categorical(wb)
xb[m,t] = xp[i, t]
end
end
return xb
end
function resample_systematic(w,M)
N = size(w,1)
bins = cumsum(exp(w)) # devec
s = collect((rand()/M+0):1/M:bins[end])
j = Array(Int64,M)
bo = 1
for i = 1:M
@inbounds for b = bo:N
if s[i] < bins[b]
j[i] = b
bo = b
break
end
end
end
return j
end
function resample_systematic!(j,w,M)
N = size(w,1)
bins = cumsum(exp(w)) # devec
s = collect((rand()/M+0):1/M:bins[end])
bo = 1
for i = 1:M
@inbounds for b = bo:N
if s[i] < bins[b]
j[i] = b
bo = b
break
end
end
end
return nothing
end
function resample_systematic_exp(w,M)
N = size(w,1)
bins = cumsum(w) # devec
s = collect((rand()/M+0):1/M:bins[end])
j = Array(Int64,M)
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.
# """
function draw_one_categorical(w)
bins = cumsum(w) # devec
s = rand()*bins[end]
for b = 1:length(bins)
if s < bins[b]
return b
end
end
end
function weighted_mean(xp,w)
N,T = size(xp)
xh = zeros(Float64,T)
for j = 1:T
@inbounds for i = 1:N
xh[j] += xp[i,j]*exp(w[i,j])
end
end
return xh
end
function weighted_mean!(xh,xp,w)
N,T = size(xp)
for j = 1:T
@inbounds for i = 1:N
xh[j] += xp[i,j]*exp(w[i,j])
end
end
return xh
end
function mode(xp,w)
N,T = size(xp)
xh = zeros(Float64,T)
for j = 1:T
i = indmax(slice(w,:,j))
xh[j] = xp[i,j]
end
return xh
end
function mode!(xh,xp,w)
N,T = size(xp)
for j = 1:T
i = indmax(slice(w,:,j))
xh[j] = xp[i,j]
end
return xh
end
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