Select Git revision
armax.jl 2.28 KiB
"""`ar(y, na; λ = 0, doplot=false)`"""
function ar(y::Vector{Float64}, na; λ = 0)
y_train, A = getARregressor(y,na)
n_points = size(A,1)
if λ == 0
w = A\y_train
else
w = (A'A + λ*eye(size(A,2)))\A'y_train
end
prediction = A*w
error = y_train - prediction
model = AR(w,na)
result = FitResult(y_train,prediction,na, λ>0?(:LS_reg) :(:LS),λ)
return model, result
end
ar(iddata::IdData, na; λ = 0, doplot=false) = ar(iddata.y, na; λ = 0)
"""arx(y, u, na, nb; λ = 0, doplot=false)"""
function arx(y::Vector{Float64}, u::VecOrMat{Float64}, na, nb; λ = 0, doplot=false)
y_train, A = getARXregressor(y,u, na, nb)
n_points = size(A,1)
if λ == 0
w = A\y_train
else
w = (A'A + λ*eye(size(A,2)))\A'y_train
end
prediction = A*w
error = y_train - prediction
si = na+1
b = Polynom[w[si:si+nb[1]-1]]
si += nb[1]
for i = 2:length(nb)
push!(b,w[si:si+nb[i]-1])
si += nb[i]
end
model = ARX(w[1:na],b,na,nb)
result = FitResult(y_train,prediction,na, λ>0?(:LS_reg) :(:LS), λ)
return model, result
end
arx(iddata::IdData, na, nb; λ = 0) = arx(iddata.y, iddata.u, na, nb; λ = 0)
function getARregressor(y::Vector{Float64},na)
A = toeplitz(y[na+1:end],y[na+1:-1:1])
y = copy(A[:,1])
A = A[:,2:end]
n_points = size(A,1)
return y,A
end
getARregressor(iddata::IdData, na) = getARXregressor(iddata.y, na)
function getARXregressor(y::Vector{Float64},u::VecOrMat{Float64}, na, nb)
assert(length(nb) == size(u,2))
m = max(na+1,maximum(nb))
n = length(y) - m+1
offs = m-na-1
A = toeplitz(y[offs+na+1:n+na+offs],y[offs+na+1:-1:1])
y = copy(A[:,1])
A = A[:,2:end]
for i = 1:length(nb)
offs = m-nb[i]
A = [A toeplitz(u[nb[i]+offs:n+nb[i]+offs-1,i],u[nb[i]+offs:-1:1+offs,i])]
end
return y,A
end
getARXregressor(iddata::IdData, na, nb) = getARXregressor(iddata.y,iddata.u, na, nb)
"""Plots the RMSE and AIC for model orders up to `n`. Useful for model selection"""
function find_na(y::Vector{Float64},n::Int)
error = zeros(n,2)
for i = 1:n
w,e = ar(y,i)
error[i,1] = rms(e)
error[i,2] = aic(e,i)
print(i,", ")
end
println("Done")
plotsub(error,"-o")
show()
end