Skip to content
Snippets Groups Projects
Commit fda9a634 authored by Fredrik Bagge Carlson's avatar Fredrik Bagge Carlson
Browse files

merge dev

parents 9810392d d411efe5
Branches
No related tags found
No related merge requests found
Showing
with 1513 additions and 204 deletions
SCS
Convex
DSP
......@@ -6,3 +6,17 @@ function PCA(W)
score = score*diagm(latent)
C,score,latent,W0
end
using MLKernels
function kernelPCA(X; α=1.0)
κ = GaussianKernel(α)
K = kernelmatrix(κ,X)
N = size(K)[1]
In = fill(1/N,(N,N))
K = K-In*K - K*In + In*K*In # Make sure mean is zero
(D,V) = eig(K)
Kpc = K*V
Kpc,D,V
end
module SystemIdentification
if !isdefined(:DEBUG); DEBUG = false; end
export
Model,LinearModel,NonLinearModel,
Network,
Polynom,PolynomMatrix,
TFdata,
AR,
ARX,
RBFARX,
......@@ -12,11 +13,12 @@ FitResult,
IdData,
# Functions
ar,arx,getARregressor,getARXregressor,find_na,
toeplitz, kalman
toeplitz, kalman, kalman_smoother, forward_kalman, PCA, plotmodel
## Fit Methods =================
:LS
:LS_reg
:L1
:LM
## Types =======================
......@@ -24,8 +26,8 @@ abstract Model
abstract LinearModel <: Model
abstract NonLinearModel <: Model
abstract Network <: NonLinearModel
typealias Polynom{T<:Real} Union(Array{T,1} , T)
typealias PolynomMatrix{T} Union(Array{Polynom{T},1},Polynom{T}, T)
typealias Polynom{T<:Real} Union{Array{T,1} , T}
typealias PolynomMatrix{T} Union{Array{Polynom{T},1},Polynom{T}, T}
......@@ -58,6 +60,11 @@ type ARX <: LinearModel
nb::Polynom{Int}
end
type TFdata <: LinearModel
P
F
end
type RBFARX <: Network
na::Int
n_centers::Int
......@@ -99,10 +106,16 @@ sse(x) = sum(x.^2)
fit(y,yh) = 100 * (1-rms(y-yh)./rms(y-mean(y)));
aic(x,d) = log(sse(x)) + 2d/length(x)
##
Base.show(fit::FitStatistics) = println("Fit RMS:$(fit.RMS), FIT:$(fit.FIT), AIC:$(fit.AIC)")
include("utilities")
include("transfer_functions")
include("idinput")
include("armax.jl")
include("kalman.jl")
include("PCA.jl")
include("kernelPCA.jl")
include("toeplitz.jl")
include("cuckooSearch.jl")
end
......@@ -7,19 +7,22 @@ function ar(y::Vector{Float64}, na; λ = 0, normtype=2, verbose=false)
if normtype == 2
if λ == 0
w = A\y_train
method = :LS
else
w = (A'A + λ*eye(size(A,2)))\A'y_train
method = :LS_reg
end
elseif normtype == 1
w = Variable(size(A,2),1)
problem = minimize(sum(abs(A*w-y_train )) + λ*norm(w))
solve!(problem, SCSSolver(verbose=Int(verbose)))
w = w.value[:]
method = :L1
end
prediction = A*w
error = y_train - prediction
model = AR(w,na)
result = FitResult(y_train,prediction,na, λ>0?(:LS_reg) :(:LS),λ)
result = FitResult(y_train,prediction,na, method,λ)
return model, result
end
ar(iddata::IdData, na; λ = 0) = ar(iddata.y, na; λ = 0)
......@@ -31,14 +34,17 @@ function arx(y::Vector{Float64}, u::VecOrMat{Float64}, na, nb; λ = 0, normtype=
if normtype == 2
if λ == 0
w = A\y_train
method = :LS
else
w = (A'A + λ*eye(size(A,2)))\A'y_train
method = :LS_reg
end
elseif normtype == 1
w = Variable(size(A,2),1)
problem = minimize(sum(abs(A*w-y_train )) + λ*norm(w))
solve!(problem, SCSSolver(verbose=Int(verbose)))
w = w.value[:]
method = :L1
end
prediction = A*w
error = y_train - prediction
......@@ -51,7 +57,7 @@ function arx(y::Vector{Float64}, u::VecOrMat{Float64}, na, nb; λ = 0, normtype=
si += nb[i]
end
model = ARX(w[1:na],b,na,nb)
result = FitResult(y_train,prediction,na, λ>0?(:LS_reg) :(:LS), λ)
result = FitResult(y_train,prediction,na, method, λ)
return model, result
end
arx(iddata::IdData, na, nb; λ = 0) = arx(iddata.y, iddata.u, na, nb; λ = 0)
......@@ -87,7 +93,7 @@ getARXregressor(iddata::IdData, na, nb) = getARXregressor(iddata.y,iddata.u, na,
"""Plots the RMSE and AIC for model orders up to `n`. Useful for model selection"""
function find_na(y::Vector{Float64},n::Int)
function find_na(y::Vector,n::Int)
error = zeros(n,2)
for i = 1:n
w,e = ar(y,i)
......@@ -96,6 +102,21 @@ function find_na(y::Vector{Float64},n::Int)
print(i,", ")
end
println("Done")
plotsub(error,"-o")
show()
scatter(error, show=true)
end
"""
plotmodel(y,m::AR)
Plots a signal `y` and the output of the model `m`
"""
function plotmodel(y,m::AR)
na = length(m.a)
y,A = getARregressor(y,na)
yh = A*m.a
error = y-yh
plot(y,c=:black)
plot(yh,c=:b)
plot(error,c=:r, title="Fitresult, AR, n_a: $na, RMSE = $(rms(error)) Fit = $(fit(y,yh))")
end
......@@ -72,6 +72,7 @@ function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=2
## Post-optimization processing
## Display all the nests
println("Total number of iterations=",N_iter);
println("f(bestnest) = $(fmin)")
squeeze(bestnest',2),fmin
end
......@@ -140,28 +141,20 @@ function empty_nests(nest,Lb,Ub,pa)
## New solution by biased/selective random walks
stepsize=rand()*(nest[randperm(n),:]-nest[randperm(n),:]);
new_nest=nest+stepsize.*K;
for j = 1:size(nest,1)
new_nest[j,:]=simplebounds(new_nest[j,:],Lb,Ub);
end
return new_nest
end
# Application of simple constraints
function simplebounds(s,Lb,Ub)
# Apply the lower bound
ns_tmp=s;
I=ns_tmp.<Lb;
ns_tmp[I]=Lb[I];
I = s.<Lb;
s[I] = Lb[I];
# Apply the upper bounds
J=ns_tmp.>Ub;
ns_tmp[J]=Ub[J];
# Update this new move
s=ns_tmp;
J = s.>Ub;
s[J] = Ub[J];
return s
end
# # ## You can replace the following by your own functions
# # # A d-dimensional objective function
# # function fobj(u)
# # ## d-dimensional sphere function sum_j=1^d (u_j-1)^2.
# # # with a minimum at (1,1, ...., 1);
# # z=sum((u-1).^2);
# # end
# # dims = 10
# # cuckoo_search(fobj,zeros(dims),-10*ones(dims),10*ones(dims))
using Devectorize
"""
`cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e5, timeout = Inf)`\n
<<<<<<< HEAD
=======
>>>>>>> dev
`n` = Number of nests (or different solutions)
`pa=0.25` Discovery rate of alien eggs/solutions
Change this if you want to get better results
Based on implementation by
@inproceedings{yang2009cuckoo,
title={Cuckoo search via L{\'e}vy flights},
author={Yang, Xin-She and Deb, Suash},
booktitle={Nature \& Biologically Inspired Computing, 2009. NaBIC 2009. World Congress on},
pages={210--214},
year={2009},
organization={IEEE}
}
http://www.mathworks.com/matlabcentral/fileexchange/29809-cuckoo-search--cs--algorithm
"""
<<<<<<< HEAD
function cuckoo_search(f,X0;Lb=-convert(Float64,Inf),Ub=convert(Float64,Inf),n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e5, timeout = Inf)
=======
function cuckoo_search(f,X0, Lb,Ub;n=25,pa=0.25, Tol=1.0e-5, max_iter = 1e3, timeout = Inf)
X00 = deepcopy(X0)
>>>>>>> dev
nd=size(X0,1);
X0t = X0'
Lb = Lb'
Ub = Ub'
if !all(isfinite(Lb))
Lb=X0t-0.99999*abs(X0t);
end
if !all(isfinite(Ub))
Ub=X0t+0.99999*abs(X0t);
end
# Random initial solutions
nest = zeros(n,nd)
<<<<<<< HEAD
nest[1,:] = X0
for i=2:n
nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb));
end
# Get the current best
fitness=10^20*ones(n,1);
fmin,bestnest,nest,fitness=get_best_nest(f,nest,nest,fitness);
=======
nest[1,:] = X0t
for i=2:n
nest[i,:]=Lb+(Ub-Lb).*rand(size(Lb));
DEBUG && @assert !any(nest[i,:] .> Ub)
DEBUG && @assert !any(nest[i,:] .< Lb)
end
# Get the current best
fitness=10^20*ones(n,1);
fmin,bestnest,nest,fitness=get_best_nest(f,nest,nest,fitness);
DEBUG && println("f(X0) = $(f(X00)), f(bestnest) = $(fmin)")
DEBUG && @assert X00 == X0
>>>>>>> dev
N_iter=0;
t0 = time()
## Starting iterations
while fmin>Tol && N_iter < max_iter
# Generate new solutions (but keep the current best)
new_nest=get_cuckoos(nest,bestnest,Lb,Ub);
fnew,best,nest,fitness=get_best_nest(f,nest,new_nest,fitness);
# Update the counter
N_iter += n;
if fnew<fmin
fmin=fnew;
bestnest=best;
end
if time()-t0 > timeout
display("Cuckoo search: timeout $(timeout)s reached ($(time()-t0)s)")
break
end
# Discovery and randomization
new_nest=empty_nests(nest,Lb,Ub,pa) ;
# Evaluate this set of solutions
fnew,best,nest,fitness=get_best_nest(f,nest,new_nest,fitness);
# Update the counter again
N_iter += n;
# Find the best objective so far
if fnew<fmin
fmin=fnew;
bestnest=best;
end
if time()-t0 > timeout
display("Cuckoo search: timeout $(timeout)s reached ($(time()-t0)s)")
break
end
end ## End of iterations
## Post-optimization processing
## Display all the nests
println("Total number of iterations=",N_iter);
<<<<<<< HEAD
=======
println("f(bestnest) = $(fmin)")
>>>>>>> dev
squeeze(bestnest',2),fmin
end
## --------------- All subfunctions are list below ------------------
## Get cuckoos by ramdom walk
function get_cuckoos(nest,best,Lb,Ub)
# Levy flights
n=size(nest,1);
# Levy exponent and coefficient
# For details, see equation (2.21), Page 16 (chapter 2) of the book
# X. S. Yang, Nature-Inspired Metaheuristic Algorithms, 2nd Edition, Luniver Press, (2010).
beta=3/2;
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
for j=1:n
s=nest[j,:];
# This is a simple way of implementing Levy flights
# For standard random walks, use step=1;
## Levy flights by Mantegna’s algorithm
u=randn(size(s))*sigma;
v=randn(size(s));
betai = 1/beta
@devec step=u./abs(v).^betai;
# In the next equation, the difference factor (s-best) means that
# when the solution is the best solution, it remains unchanged.
stepsize=0.01*step.*(s-best);
# Here the factor 0.01 comes from the fact that L/100 should the typical
# step size of walks/flights where L is the typical lenghtscale;
# otherwise, Levy flights may become too aggresive/efficient,
# which makes new solutions (even) jump out side of the design domain
# (and thus wasting evaluations).
# Now the actual random walks or flights
s=s+stepsize.*randn(size(s));
# Apply simple bounds/limits
nest[j,:]=simplebounds(s,Lb,Ub);
<<<<<<< HEAD
end
nest
=======
DEBUG && @assert !any(nest[j,:] .> Ub)
DEBUG && @assert !any(nest[j,:] .< Lb)
end
return nest
>>>>>>> dev
end
## Find the current best nest
function get_best_nest(f,nest,newnest,fitness)
# Evaluating all new solutions
for j=1:size(nest,1)
fnew=f(squeeze(newnest[j,:]',2));
if fnew<=fitness[j]
fitness[j]=fnew;
nest[j,:]=newnest[j,:];
end
end
# Find the current best
(fmin,K) = findmin(fitness) ;
best=nest[K,:];
<<<<<<< HEAD
fmin,best,nest,fitness
=======
return fmin,best,nest,fitness
>>>>>>> dev
end
## Replace some nests by constructing new solutions/nests
function empty_nests(nest,Lb,Ub,pa)
# A fraction of worse nests are discovered with a probability pa
n=size(nest,1);
# Discovered or not -- a status vector
K=rand(size(nest)).>pa;
# In the real world, if a cuckoo’s egg is very similar to a host’s eggs, then
# this cuckoo’s egg is less likely to be discovered, thus the fitness should
# be related to the difference in solutions. Therefore, it is a good idea
# to do a random walk in a biased way with some random step sizes.
## New solution by biased/selective random walks
stepsize=rand()*(nest[randperm(n),:]-nest[randperm(n),:]);
new_nest=nest+stepsize.*K;
<<<<<<< HEAD
=======
for j = 1:size(nest,1)
new_nest[j,:]=simplebounds(new_nest[j,:],Lb,Ub);
end
return new_nest
>>>>>>> dev
end
# Application of simple constraints
function simplebounds(s,Lb,Ub)
# Apply the lower bound
<<<<<<< HEAD
ns_tmp=s;
I=ns_tmp.<Lb;
ns_tmp[I]=Lb[I];
# Apply the upper bounds
J=ns_tmp.>Ub;
ns_tmp[J]=Ub[J];
# Update this new move
s=ns_tmp;
=======
I = s.<Lb;
s[I] = Lb[I];
# Apply the upper bounds
J = s.>Ub;
s[J] = Ub[J];
return s
>>>>>>> dev
end
# # ## You can replace the following by your own functions
# # # A d-dimensional objective function
# # function fobj(u)
# # ## d-dimensional sphere function sum_j=1^d (u_j-1)^2.
# # # with a minimum at (1,1, ...., 1);
# # z=sum((u-1).^2);
# # end
# # dims = 10
# # cuckoo_search(fobj,zeros(dims),-10*ones(dims),10*ones(dims))
"""
idinput(N, class = "prbs"; band = 1)
Create a input signal for a sys.id. experiment, currently only supports prbs signals.
"""
function idinput(N, class = "prbs"; band = 1)
output = zeros(N)
if lowercase(class) == "prbs"
......@@ -18,4 +23,3 @@ function idinput(N, class = "prbs"; band = 1)
return output
end
function kalman(R1,R2,theta, y, A, P)
# This is really a Kalman filter, not RLS
ATP = A'*P;
K = (P*A)/(R2+ATP*A);
P = P - (P*A*ATP)./(R2 + ATP*A) + R1;
yp = A'*theta
e = (y-yp)[1];
red = 1;
# if abs(e) > 0.025
# red = 0.2;
# end
theta = theta + K*e*red;
theta, P, e, yp[1]
"""
`kalman(R1,R2,theta, y, A, P)`
One dimensional Kalman filter for parameter estimates
"""
function kalman1(R1,R2,theta, y, A, P)
ATP = A'P
K = (P*A)/(R2+ATP*A)
P = P - (P*A*ATP)./(R2 + ATP*A) + R1
yp = (A'theta)[1]
e = y-yp
theta = theta + K*e
return theta, P, e, yp
end
"""
`kalman(R1,R2,R12,theta, y, A, P)`
General Kalman filter for parameter estimates
"""
function kalman(R1,R2,R12,theta, y, A, P)
ATP = A'P
K = (P*A+R12)/(R2+ATP*A)
P = P - (P*A+R12)/(R2 + ATP*A)*(ATP+R12') + R1
yp = (A'theta)[1]
e = y-yp
theta = theta + K*e
return theta, P, e, yp
end
function forward_kalman(y,A,R1,R2, P0)
na = size(R1,1)
N = length(y)
xkk = zeros(na,N);
Pkk = zeros(na,na,N)
xkk[:,1] = A\y;
Pkk[:,:,1] = P0;
xk = xkk
Pk = Pkk
i = 1
Hk = A[i,:]
e = y[i]-Hk*xk[:,i]
S = Hk*Pk[:,:,i]*Hk' + R2
K = (Pk[:,:,i]*Hk')/S
xkk[:,i] = xk[:,i] + K*e
Pkk[:,:,i] = (I - K*Hk)*Pk[:,:,i]
for i = 2:N
Fk = 1
Hk = A[i,:]
xk[:,i] = Fk*xkk[:,i-1]
Pk[:,:,i] = Fk*Pkk[:,:,i-1]*Fk' + R1
e = y[i]-Hk*xk[:,i]
S = Hk*Pk[:,:,i]*Hk' + R2
K = (Pk[:,:,i]*Hk')/S
xkk[:,i] = xk[:,i] + K*e
Pkk[:,:,i] = (I - K*Hk)*Pk[:,:,i]
end
return xkk,xk,Pkk,Pk
end
using Debug
# """A kalman parameter smoother"""
function kalman_smoother(y, A, R1, R2, P0)
na = size(R1,1)
N = length(y)
xkk,xk,Pkk,Pk = forward_kalman(y,A,R1,R2, P0)
xkn = zeros(xkk)
Pkn = zeros(Pkk)
for i = N-1:-1:1
Fk = 1
Hk = A[i,:]'
C = Pkk[:,:,i]/Pk[:,:,i+1]
xkn[:,i] = xkk[:,i] + C*(xkn[:,i+1] - xk[:,i+1])
Pkn[:,:,i] = Pkk[:,:,i] + C*(Pkn[:,:,i+1] - Pk[:,:,i+1])*C'
end
return xkn
end
function kernelPCA(X; α=1.0)
κ = GaussianKernel(α)
K = kernelmatrix(κ,X)
N = size(K)[1]
In = ones(N,N)./N
K = K-In*K - K*In + In*K*In
(D,V) = eig(K)
Kpc = K*V
Kpc,D,V
end
"""
Calculates the Tikhonov regularized estimate min. ||Ax-b|| + λ||x|| using SVD
for a vector with different values of λ
`x, normE, normX = dsvd(A,U,S,V,b,λ)`
`x, normE, normX = dsvd(A,b,λ)`
"""
function dsvd(A,U,S,V,b,λ)
if minimum(λ) < 0
error("Illegal regularization parameter λ")
end
UTb = U'b
N = size(U,1)
n = size(V,1)
nl = length(λ)
x = zeros(n,nl)
normE = zeros(nl)
normX = zeros(nl)
print("dsvd: ")
for (i,λi) in enumerate(λ)
x[:,i] = V*diagm(S./(S.^2+λi))*UTb
normX[i] = norm(x[:,i])
normE[i] = norm(A*x[:,i] - b)
print("[$(round(i/length(λ),2)),$(round(λi,2))] ")
end
println("Done ")
return x, normE, normX
end
function dsvd(A,b,λ)
if minimum(λ) < 0
error("Illegal regularization parameter λ")
end
U,S,V = svd(A)
dsvd(A,U,S,V,b,λ)
end
"""
Lcurve(normE, normX, λ)
Plots the L-curve. normE and normX are obtained from, e.g., `dsvd`
"""
function Lcurve(normE, normX, λ)
plot(normE,normX,xscale=:log10,yscale=:log10,m=:o, xlabel="RMSE", ylabel="||k||", title="L-curve")
annotations = [(normE[i],normX[i],"λ=$(round(λ[i],8))") for i in 1:length(λ)]
annotate!(annotations)
end
# module Tmp
using Devectorize
include("pf_bootstrap.jl")
include("pf_smoother_FFBSi.jl")
include("pf_PMMH.jl")
include("pf_PG.jl")
include("resample.jl")
include("utilities.jl")
const σw0 = 2.0
const σw = 0.32
const σv = 1.0
const theta0 = [0.5, 25, 8]
s2piσv = log(sqrt(2*pi) * σv)
function f(x::Vector, 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
end
end
f(x::Float64, t::Int64) = 0.5*x + 25*x/(1+x^2) + 8*cos(1.2*(t-1))
f(xn,x,t::Int64) = begin
c = 8*cos(1.2*(t-1))
@inbounds for i = 1:length(x)
xn[i] = 0.5*x[i] + 25*x[i]./(1+x[i]^2) + c + σw*randn()
end
end
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()
end
return x
end
f_sample(x::Float64, t::Int64) = 0.5*x + 25*x/(1+x^2) + 8*cos(1.2*(t-1)) + σw*randn()
function f_density(xb, xp, t)
x = f(xp,t)-xb
exp(-0.5 * (x/σw)^2)
end
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)
end
end
## Parameter estimation
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]*x1[i] + theta[2]*x1[i]./(1+x1[i]^2) + c + σw*randn()
end
return x
end
f_theta_sample(x::Float64, t, theta) = theta[1]*x + theta[2]*x./(1+x^2) + theta[3]*cos(1.2*(t-1)) + σw*randn()
function f_density_theta(xb, xp, t, theta)
x = f_theta(xp,t, theta)-xb
exp(-0.5 * (x/σw)^2)
end
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 - s2piσv
end
return w
end
draw_theta_random(theta) = theta + 0.01*randn(size(theta0)).*abs(theta0)
draw_theta0() = theta0+ 0.2*randn(size(theta0)).*abs(theta0)
function draw_theta_PG(theta, x)
t = 0:length(x)-1
[x x./(1+x.^2) cos(1.2*t)][1:end-1,:]\x[2:end]# + draw_theta_random(theta)
end
rms(x) = sqrt(mean(x.^2))
function test_pf()
particle_count = [5 15 50 150 500 3000 10_000]
time_steps = [100, 2000]
RMSE = zeros(length(particle_count),length(time_steps)) # Store the RMS errors
RMSE_FFBSi = zeros(length(particle_count),length(time_steps))
propagated_particles = 0
M = 10
for (Ti,T) in enumerate(time_steps)
xh = Array(Float64,T)
for (Ni, N) in enumerate(particle_count)
montecarlo_runs = maximum(particle_count)*maximum(time_steps) / T / N*10
# montecarlo_runs = 1
xp = Array(Float64,(N,T))
w = Array(Float64,(N,T))
x = Array(Float64,T)
y = Array(Float64,T)
for mc_iter = 1:montecarlo_runs
x[1] = σw*randn()
y[1] = σv*randn()
for t = 1:T-1
x[t+1] = f_sample(x[t],t)
y[t+1] = 0.05x[t+1]^2 + σv*randn()
end # t
pf_bootstrap!(xp,w, y, N, g_density, f_sample, σw0 )
xb = pf_smoother_FFBSi( w,xp, M, f_density )
# weighted_mean!(xh,xp,w)
mode!(xh,xp,w)
RMSE[Ni,Ti] += rms(x-xh)
RMSE_FFBSi[Ni,Ti] += rms(x-mean(xb,1)')
end # MC
RMSE[Ni,Ti] /= montecarlo_runs
RMSE_FFBSi[Ni,Ti] /= montecarlo_runs
propagated_particles += montecarlo_runs*N*T
@show N
end # N
@show T
end # T
println("Propagated $propagated_particles particles")
#
return RMSE, RMSE_FFBSi
end
function test_PMMH()
particle_count = [100]
time_steps = [200]
R = 30_000
for (Ti,T) in enumerate(time_steps)
for (Ni, N) in enumerate(particle_count)
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))
@sync @parallel for mc_iter = 1:montecarlo_runs
x[1] = σw*randn()
y[1] = σv*randn()
for t = 1:T-1
x[t+1] = f_sample(x[t],t)
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, 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
end # T
println("Done test_PMMH")
# return theta
end
function test_PG()
particle_count = [10]
time_steps = [500]
R = 10_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
x = Array(Float64,T)
y = Array(Float64,T)
theta_PG = SharedArray(Float64,(size(theta0,1),R,montecarlo_runs))
@sync @parallel for mc_iter = 1:montecarlo_runs
x[1] = σw*randn()
y[1] = σv*randn()
for t = 1:T-1
x[t+1] = f_sample(x[t],t)
y[t+1] = 0.05x[t+1]^2 + σv*randn()
end # t
theta0i = draw_theta0()
theta_PG[:,:,mc_iter] = pf_PG(y, N, R, pf_smoother_FFBSi!, f_theta_sample, draw_theta_PG, f_density_theta, g_density_theta, σw0, theta0i )
end # MC
newplot(theta_PG[:,:,1]'); hold(true)
for mc_iter = 2:montecarlo_runs
plot(theta_PG[:,:,mc_iter]')
end
hold(false); title("PG")
th = Array(theta_PG)
figure(); pp = plothist(th[1,:,:][:],2000); hold(true)
plothist(th[2,:,:][:],2000)
plothist(th[3,:,:][:],2000)
hold(false); display(pp)
return theta_PG
@show N
end # N
@show T
end # T
println("Done test_PMMH")
# return theta
end
function test_PMMH_PG()
particle_count = [100]
time_steps = [200]
R = 5_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
x = Array(Float64,T)
y = Array(Float64,T)
theta_PMMH = SharedArray(Float64,(size(theta0,1),R,montecarlo_runs))
theta_PG = SharedArray(Float64,(size(theta0,1),R,montecarlo_runs))
@sync @parallel for mc_iter = 1:montecarlo_runs
x[1] = σw*randn()
y[1] = σv*randn()
for t = 1:T-1
x[t+1] = f(x[t],t) + σw*randn()
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, draw_theta_random, f_density_theta, g_density_theta, σw0, theta0i)
theta_PG[:,:,mc_iter] = pf_PG(y, N, R, pf_smoother_FFBSi!, f_theta, draw_theta_PG, f_density_theta, g_density_theta, σw0, 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")
newplot(theta_PG[:,:,1]'); hold(true)
for mc_iter = 2:montecarlo_runs
plot(theta_PG[:,:,mc_iter]')
end
hold(false); title("PG")
@show N
end # N
@show T
end # T
println("Done test_PMMH")
# return theta
end
function plotting(RMSE)
particle_count = [5 15 50 150 500 3000 10_000]
time_steps = [100, 2000]
nT = length(time_steps)
p = loglog(particle_count,RMSE,"o")
legend_strings = Array(String,nT)
for i = 1:nT
legend_strings[i] = "$(time_steps[i]) time steps"
end
title("RMS errors vs Number of particles")
@windows_only display(p)
@linux_only PYPLOT && legend(legend_strings)
# ProfileView.view()
end
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)
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_PG(y, N, R, pf_smoother, 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))
xb = Array(Float64,1,T)
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_density(xb,xp,t) = f_density_theta(xb, xp, t, theta0)
f(x,t) = f_theta(x,t,theta0)
pf_bootstrap!(xp, w, y, N, g_density, f)
pf_smoother(xb, w,xp, 1, f_density )
breakPoints = 10
# noiseMultiplier = 0.1^(1/R)
Nstart = N
breakTime = floor(Int64, R / breakPoints)
particleIncrease = 10*N
thetaP = draw_theta(theta0,xb')
f(x,t) = f_theta(x,t,thetaP)
g_density(y_t, x_t,w) = g_density_theta(y_t, x_t,w, thetaP)
f_density(xb,xp,t) = f_density_theta(xb, xp, t, thetaP)
for r = 2:R
try
pf_CSMC!(xp, w, y, N, g_density, f, xb)
pf_smoother(xb, w,xp, 1, f_density )
# @show exp(Zp-Z[r-1])
thetaP = draw_theta(theta[:,r-1],xb')
theta[:,r] = thetaP
catch ex
display(ex)
@show r
@show thetaP
error("Stopping")
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
function pf_PMMH(y, N, R, pf, f_theta, draw_theta, f_density_theta, g_density_theta, 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,x1,t) = f_theta(x,x1,t,theta0)
pf(xp, w, wnn, y, N, g_density, f)
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,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
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)
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
# module Tmp
using Devectorize
using StatsBase
function pf_bootstrap!(xp,w, y, N, g_density, f)
T = length(y)
xp[:,1] = 2*σw*randn(N)
wT = slice(w,:,1)
xT = slice(xp,:,1)
fill!(wT,-log(N))
g_density(y[1],xT, wT)
wT -= logsumexp(wT)
j = Array(Int64,N)
# tw = Float64(N)
for t = 2:T
xT1 = xT
xT = slice(xp,:,t)
wT1 = wT
wT = slice(w,:,t)
# Resample
if t%2 == 0
resample_systematic!(j,wT1,N)
f(xT,xT1[j],t-1)
fill!(wT,-log(N))
else # Resample not needed
f(xT,xT1,t-1)
copy!(wT,wT1)
end
# Time update
g_density(y[t],xT, wT)
# Normalize weights
wT -= logsumexp(wT)
# tw = 1/sum(w(:,t).*2)
end
return nothing
end
function pf_bootstrap(y, N, g_density, f)
T = length(y)
xp = Array(Float64,(N,T))
w = Array(Float64,(N,T))
xp[:,1] = 2*σw*randn(N)
wT = slice(w,:,1)
xT = slice(xp,:,1)
fill!(wT,-log(N))
g_density(y[1],xT, wT)
wT -= logsumexp(wT)
j = Array(Int64,N)
# tw = Float64(N)
for t = 2:T
xT1 = xT
xT = slice(xp,:,t)
wT1 = wT
wT = slice(w,:,t)
# Resample
if t%2 == 0
resample_systematic!(j,wT1,N)
f(xT,xT1[j],t-1)
fill!(wT,-log(N))
else # Resample not needed
f(xT,xT1,t-1)
copy!(wT,wT1)
end
# Time update
g_density(y[t],xT, wT)
# Normalize weights
wT -= logsumexp(wT)
# tw = 1/sum(w(:,t).*2)
end
return xp,w
end
function pf_bootstrap_nn(xp, w, wnn, y, N, g_density, f)
T = length(y)
xp[:,1] = 2*σw*randn(N)
wT = slice(w,:,1)
xT = slice(xp,:,1)
fill!(wT,-log(N))
g_density(y[1],xT, wT)
wnn[:,1] = wT
wT -= logsumexp(wT)
j = Array(Int64,N)
# tw = Float64(N)
for t = 2:T
xT1 = xT
xT = slice(xp,:,t)
wT1 = wT
wT = slice(w,:,t)
# Resample
if t%2 == 0
resample_systematic!(j,wT1,N)
f(xT,xT1[j],t-1)
fill!(wT,-log(N))
else # Resample not needed
f(xT,xT1,t-1)
copy!(wT,wT1)
end
# Time update
g_density(y[t],xp[:,t], wT)
wnn[:,t] = wT
# Normalize weights
#some strange error with this normalization method, I totally have no clue why it's working elsewhere but not here
wT -= logsumexp(wT)
# @assert isapprox(nc,log(sum(exp(wnn[:,t]))))
# tw = 1/sum(w(:,t).*2)
end
return xp,w, wnn
end
function pf_aux_nn(xp, w, wnn, y, N, g_density, f)
T = length(y)
xp[:,1] = 2*σw*randn(N)
w0 = fill(log(1/N),N)
xp[:,1] = f(zeros(N),1)
wnn[:,1] = w0 + g_density(y[1],xp[:,1])
@devec w[:,1] = wnn[:,1] - log(sum(exp(wnn[:,1])))
j = Array(Int64,N)
# tw = Float64(N)
for t = 2:T
xpT = xp[:,t-1]
f(xpT,t-1) # This is now incorrect. particles should be propagated without noise in this step.
wT = copy(w[:,t-1])
lambda = g_density(y[t],xpT)
wT += lambda
wTe = exp(wT)
wTe /= sum(wTe)
# Resample
if t%2 == 0
resample_systematic_exp!(j,wTe)
xpT = xp[j,t-1]
wT = copy(w0)
lambdaT = lambda[j]
else # Resample not needed
xpT = xp[:,t-1]
# wT = copy(w[:,t-1])
lambdaT = lambda
end
# Time update
f(xpT,t-1)
xp[:,t] = xpT
wT += g_density(y[t],xp[:,t])
wnn[:,t] = wT - lambdaT
# Normalize weights
#some strange error with this normalization method, I totally have no clue why it's working elsewhere but not here
# # offset = maximum(wT)
# # normConstant = 0.0
# # for i = 1:N
# # normConstant += exp(wT[i]-offset)
# # end
# # w[:,t] = wT - log(normConstant)+offset
w[:,t] = wnn[:,t] - log(sum(exp(wnn[:,t])))
# tw = 1/sum(w(:,t).*2)
end
return xp,w, wnn
end
function pf_CSMC!(xp,w, y, N, g_density, f, xc )
xp[:,1] = 2*σw*randn(N)
xp[N,1] = xc[1]
T = length(y)
wT = slice(w,:,1)
xT = slice(xp,:,1)
fill!(wT,-log(N))
g_density(y[1],xT, wT)
wT -= logsumexp(wT)
j = Array(Int64,N)
# tw = Float64(N)
for t = 2:T
xT1 = xT
xT = slice(xp,:,t)
wT1 = wT
wT = slice(w,:,t)
# Resample
if t%2 == 0
resample_systematic!(j,wT1,N)
f(xT,xT1[j],t-1)
fill!(wT,-log(N))
else # Resample not needed
f(xT,xT1,t-1)
copy!(wT,wT1)
end
xp[N,t] = xc[t]
# Time update
g_density(y[t],xT, wT)
# Normalize weights
wT -= logsumexp(wT)
# tw = 1/sum(w(:,t).*2)
end
return nothing
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 pf_smoother_FFBSi!( xb, w,xp, M, f_density )
(N,T) = size(w)
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
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
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.
# """
function draw_one_categorical(w)
bins = cumsum(w) # devec
s = rand()*bins[end]
midpoint = round(Int64,length(bins)/2)
if s < bins[midpoint]
for b = 1:midpoint
if s < bins[b]
return b
end
end
else
for b = midpoint:length(bins)
if s < bins[b]
return b
end
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
function getΨVec!(Ψ,σvec,phi,centers,w=0;period = 0, normalized=false)
if period > 0
psi(Phi,c,σ) = exp((cos((Phi-c)*period/(2*pi))-1)/σ);
else
function psi!(ψ,j,Phi,c,σ)
@simd for i = 1:size(ψ,1)
ψ[i,j] = exp(-(c-Phi[i])^2*σ)[1]
end
end
end
if w == 0
w = zeros(size(Ψ,2))
end
σvec = σvec[:]
centers = centers[:]
Nsignals = 1
(Nbasis, ) = size(centers);
= zeros(size(Ψ,1),size(centers,1)*2)
if normalized
for i = 1:Nbasis
psi!(Ψ,i,phi,centers[i],σvec[i]);
@devec Ψ[:,i] = Ψ[:,i]./sum(Ψ[:,i])
[:,i] = -Ψ[:,i].*2.*(centers[i] - phi).*σvec[i]
[:,i+Nbasis] = -Ψ[:,i].*(centers[i] - phi).^2
end
else
for i = 1:Nbasis
psi!(Ψ,i, phi,centers[i],σvec[i]);
[:,i] = -Ψ[:,i].*2.*(centers[i] - phi).*σvec[i]
[:,i+Nbasis] = -Ψ[:,i].*(centers[i] - phi).^2
end
end
return
# if normalized
# Ψ = Ψ./repmat(sum(Ψ,2),1,size(Ψ,2))
# end
end
function getΨVecNd!(Ψ, phi , params, w=0; period = 0, normalized=true)
Nbasis = size(Ψ,2)
Nsignals = size(phi,2)
if period > 0
psi(Phi,c,σ) = exp((cos((c-Phi)*period/(2*pi))-1)/σ);
else
if Nsignals >= 1
function psi(Phi,c,σ)
# Denna jävla funktionen allokerar en massa minne!
ret = exp(- sum(((c-Phi).^2).*σ) )
end
else
function psi(Phi,c,σ)
# Denna jävla funktionen allokerar en massa minne!
ret = exp(- ((Phi-c)^2*(1/σ))[1,1] )
end
end
end
if w == 0
w = zeros(Nbasis)
end
# figure(); pp= imagesc(Centers); display(pp);
# figure(); pp= imagesc(Σ); display(pp);
Σoffset = convert(Int64,size(params,1)/2);
N = size(phi,1)
= zeros(N,size(params,1))
if normalized
s = 1e-10*ones(size(phi,1))
i1 = 1
for i = 1:Nbasis
for j = 1:size(phi,1)
Ψ[j,i] = psi(phi[j,:]',params[i1:(i1+Nsignals-1)],params[(Σoffset+i1):(Σoffset+i1+Nsignals-1)])
s[j] += Ψ[j,i]
end
i1 += Nsignals
end
for j = 1:size(phi,1)
@devec Ψ[j,:] = Ψ[j,:]./s[j]
end
i1 = 1
for i = 1:Nbasis
for j = 1:N
c = params[i1:(i1+Nsignals-1)]
s = params[(Σoffset+i1):(Σoffset+i1+Nsignals-1)]
x = phi[j,:]'
[j,i1:(i1+Nsignals-1)] = (Ψ[j,i]*w[i]).*(2*(x-c).*s)
[j,(Σoffset+i1):(Σoffset+i1+Nsignals-1)] = -(Ψ[j,i]*w[i]).*((c-x).^2 )
end
i1 += Nsignals
end
else
i1 = 1
for i = 1:Nbasis
for j = 1:N
c = params[i1:(i1+Nsignals-1)]
s = params[(Σoffset+i1):(Σoffset+i1+Nsignals-1)]
x = phi[j,:]'
Ψ[j,i] = psi(x,c,s)[1]
[j,i1:(i1+Nsignals-1)] = (Ψ[j,i]*w[i]).*(2*(x-c).*s)
[j,(Σoffset+i1):(Σoffset+i1+Nsignals-1)] = -(Ψ[j,i]*w[i]).*((c-x).^2 )
end
i1 += Nsignals
end
end
return
end
function getCenters(centers,σvec)
(nbasis,Nsignals) = size(centers)
Nbasis::Int64 = nbasis^Nsignals
Centers = zeros(Nsignals, Nbasis)
Σ = zeros(Nsignals, Nbasis)
v = Nbasis
h = 1
for i = 1:Nsignals
v = convert(Int64,v / nbasis)
Centers[i,:] = vec(repmat(centers[:,i]',v,h))'
Σ[i,:] = vec(repmat(σvec[:,i]',v,h))'
h *= nbasis
end
[vec(Centers), vec(Σ)]
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment