diff --git a/src/kalman_base.jl b/src/kalman_base.jl
index e31013f97b809c1577f40fb629ecfd7e13c2301f..eaca018cbf2d98c6acaa17e19aea0e048a1da24c 100644
--- a/src/kalman_base.jl
+++ b/src/kalman_base.jl
@@ -298,13 +298,13 @@ function get_updated_Pstartt!(Pstartt::AbstractMatrix{T}, Pstar::AbstractMatrix{
end
# v = y - Z*a -- basic
-function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{T}, z::AbstractVecOrMat{T}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
+function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractVecOrMat{T}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
copyto!(v, 1, y, iy, ny)
gemv!('N', -1.0, z, a, 1.0, v)
end
# v = y - Z*a -- basic -- univariate
-function get_v!(Y::AbstractVecOrMat{T}, Z::AbstractVecOrMat{T}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer}
+function get_v!(Y::AbstractVecOrMat{Union{T, Missing}}, Z::AbstractVecOrMat{T}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer}
v = Y[i]
@inbounds @simd for j = 1:length(a)
v -= Z[i, j]*a[j]
@@ -313,37 +313,37 @@ function get_v!(Y::AbstractVecOrMat{T}, Z::AbstractVecOrMat{T}, a::AbstractVecto
end
# v = y - a[z] -- Z selection matrix
-function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{T}, z::AbstractVector{U}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
+function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
copyto!(v, 1, y, iy, ny)
az = view(a,z)
v .= v .- az
end
# v = y - a[z] -- Z selection matrix -- univariate
-function get_v!(y::AbstractVecOrMat{T}, z::AbstractVector{U}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer}
+function get_v!(y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer}
return y[i] - a[z[i]]
end
# v = y - Z*a -- missing observations
-function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{T}, z::AbstractMatrix{T}, a::AbstractVector{T}, t::U, pattern::Vector{U}) where {T <: AbstractFloat, U <: Integer}
+function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractMatrix{T}, a::AbstractVector{T}, t::U, pattern::Vector{U}) where {T <: AbstractFloat, U <: Integer}
v .= view(y, pattern, t)
gemv!('N', -1.0, z, a, 1.0, v)
end
# v = y - a[z] -- Z selection matrix and missing variables
-function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{T}, z::AbstractVector{U}, a::AbstractVector{T}, t::U, pattern::Vector{Int64}) where {T <: AbstractFloat, U <: Integer}
+function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{Union{T, Missing}}, z::AbstractVector{U}, a::AbstractVector{T}, t::U, pattern::Vector{Int64}) where {T <: AbstractFloat, U <: Integer}
v .= view(y, pattern, t) .- view(a, z)
end
# v = y - c - Z*a -- basic
-function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{T}, c::AbstractVector{T}, z::AbstractArray{T}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
+function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractVector{T}, z::AbstractArray{T}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
copyto!(v, 1, y, iy, ny)
v .-= c
gemm!('N', 'N', -1.0, z, a, 1.0, v)
end
# v = y - c - Z*a -- basic -- univariate
-function get_v!(Y::AbstractVecOrMat{T}, c::AbstractVector{T}, Z::AbstractVecOrMat{T}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer}
+function get_v!(Y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractVector{T}, Z::AbstractVecOrMat{T}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer}
v = Y[i] - c[i]
@inbounds @simd for j = 1:length(a)
v -= Z[i, j]*a[j]
@@ -352,25 +352,25 @@ function get_v!(Y::AbstractVecOrMat{T}, c::AbstractVector{T}, Z::AbstractVecOrMa
end
# v = y - c - a[z] -- Z selection matrix
-function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{T}, c::AbstractArray{T}, z::AbstractVector{U}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
+function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractArray{T}, z::AbstractVector{U}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
copyto!(v, 1, y, iy, ny)
az = view(a,z)
v .-= c .+ az
end
# v = y - c - a[z] -- Z selection matrix -- univariate
-function get_v!(y::AbstractVecOrMat{T}, c::AbstractVector{T}, z::AbstractVector{U}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer}
+function get_v!(y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractVector{T}, z::AbstractVector{U}, a::AbstractVector{T}, i::U) where {T <: AbstractFloat, U <: Integer}
return y[i] - c[i] - a[z[i]]
end
# v = y - c - Z*a -- missing observations
-function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{T}, c::AbstractArray{T}, z::AbstractArray{T}, a::AbstractArray{T}, t::U, pattern::Vector{U}) where {T <: AbstractFloat, U <: Integer}
+function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractArray{T}, z::AbstractArray{T}, a::AbstractArray{T}, t::U, pattern::Vector{U}) where {T <: AbstractFloat, U <: Integer}
v .= view(y, pattern, t) .- view(c, pattern)
gemm!('N', 'N', -1.0, z, a, 1.0, v)
end
# v = y - c - a[z] -- Z selection matrix and missing variables
-function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{T}, c::AbstractArray{T}, z::AbstractVector{U}, a::AbstractArray{T}, t::U, pattern::Vector{Int64}) where {T <: AbstractFloat, U <: Integer}
+function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{Union{T, Missing}}, c::AbstractArray{T}, z::AbstractVector{U}, a::AbstractArray{T}, t::U, pattern::Vector{Int64}) where {T <: AbstractFloat, U <: Integer}
v .= view(y, pattern, t) .- view(c, pattern) .- view(a, z)
end
diff --git a/src/kalman_filter.jl b/src/kalman_filter.jl
index a99a66d4e4f588c26bae3b283a37f6334138862e..5e2f45296d8b87e4498c42aeb7f9aba3bed4fa0a 100644
--- a/src/kalman_filter.jl
+++ b/src/kalman_filter.jl
@@ -94,7 +94,7 @@ end
KalmanFilterWs(ny, ns, np, nobs) = KalmanFilterWs{Float64, Int64}(ny, ns, np, nobs)
-function kalman_filter!(Y::AbstractArray{U},
+function kalman_filter!(Y::AbstractArray{Union{U, Missing}},
c::AbstractArray{U},
Z::AbstractArray{W},
H::AbstractArray{U},
@@ -140,7 +140,7 @@ function kalman_filter!(Y::AbstractArray{U},
pattern = data_pattern[t]
ndata = length(pattern)
vc = changeC ? view(c, :, t) : view(c, :)
- ws.csmall .= view(vc, pattern)
+ # ws.csmall .= view(vc, pattern)
vZsmall = get_vZsmall(ws.Zsmall, ws.iZsmall, Z, pattern, ndata, ny, t)
vH = changeH ? view(H, :, :, t) : view(H, :, :)
vT = changeT ? view(T, :, :, t) : view(T, :, :)
@@ -164,32 +164,29 @@ function kalman_filter!(Y::AbstractArray{U},
vvH = view(vH, pattern, pattern)
vZP = view(ws.ZP, 1:ndata, :)
vcholF = view(ws.cholF, 1:ndata, 1:ndata, t)
+ vcholH = view(ws.cholH, 1:ndata, 1:ndata)
# v = Y[:,t] - c - Z*a
get_v!(vv, Y, vc, vZsmall, va, t, pattern)
- if !steady
- # F = Z*P*Z' + H
- get_F!(vF, vZP, vZsmall, vP, vvH)
- info = get_cholF!(vcholF, vF)
- if info != 0
- # F is near singular
- if !cholHset
- get_cholF!(ws.cholH, vvH)
- cholHset = true
- end
- ws.lik[t] = ndata*l2pi + univariate_step!(vatt, va1, vPtt, vP1, Y, t, c, ws.Zsmall, vvH, d, T, ws.QQ, va, vP, ws.kalman_tol, ws, pattern)
- t += 1
- continue
+ # F = Z*P*Z' + H
+ get_F!(vF, vZP, vZsmall, vP, vvH)
+ info = get_cholF!(vcholF, vF)
+ if info != 0
+ # F is near singular
+ if !cholHset
+ get_cholF!(vcholH, vvH)
+ cholHset = true
end
+ ws.lik[t] = ndata*l2pi + univariate_step!(vatt, va1, vPtt, vP1, Y, t, c, ws.Zsmall, vvH, d, T, ws.QQ, va, vP, ws.kalman_tol, ws, pattern)
+ t += 1
+ continue
end
# iFv = inv(F)*v
get_iFv!(viFv, vcholF, vv)
ws.lik[t] = ndata*l2pi + log(det_from_cholesky(vcholF)) + LinearAlgebra.dot(vv, viFv)
if t <= last
- if !steady
- # K = iF*Z*P
- get_K!(vK, vZP, vcholF)
- end
+ # K = iF*Z*P
+ get_K!(vK, vZP, vcholF)
# att = a + K'*v
get_updated_a!(vatt, va, vK, vv)
# a = d + T*att
@@ -289,7 +286,7 @@ end
DiffuseKalmanFilterWs(ny, ns, np, nobs) = DiffuseKalmanFilterWs{Float64, Int64}(ny, ns, np, nobs)
-function diffuse_kalman_filter_init!(Y::AbstractArray{U},
+function diffuse_kalman_filter_init!(Y::AbstractArray{Union{U, Missing}},
c::AbstractArray{U},
Z::AbstractArray{W},
H::AbstractArray{U},
@@ -311,7 +308,6 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
data_pattern::Vector{Vector{V}}) where {U <: AbstractFloat,
V <: Integer,
W <: Real}
-
changeC = ndims(c) > 1
changeH = ndims(H) > 2
changeD = ndims(d) > 1
@@ -337,7 +333,7 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
pattern = data_pattern[t]
ndata = length(pattern)
vc = changeC ? view(c, :, t) : view(c, :)
- ws.csmall .= view(vc, pattern)
+ # ws.csmall .= view(vc, pattern)
vZsmall = get_vZsmall(ws.Zsmall, ws.iZsmall, Z, pattern, ndata, ny, t)
vH = changeH ? view(H, :, :, t) : view(H, :, :)
vT = changeT ? view(T, :, :, t) : view(T, :, :)
@@ -419,7 +415,7 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
return t
end
-function diffuse_kalman_filter!(Y::AbstractArray{U},
+function diffuse_kalman_filter!(Y::AbstractArray{Union{U, Missing}},
c::AbstractArray{U},
Z::AbstractArray{W},
H::AbstractArray{U},
@@ -451,7 +447,7 @@ function diffuse_kalman_filter!(Y::AbstractArray{U},
return -0.5*sum(vlik)
end
-function diffuse_kalman_filter!(Y::AbstractArray{U},
+function diffuse_kalman_filter!(Y::AbstractArray{Union{U, Missing}},
c::AbstractArray{U},
Z::AbstractArray{W},
H::AbstractArray{U},
diff --git a/src/kalman_likelihood.jl b/src/kalman_likelihood.jl
index ed34d7cdf74fcdaf1e80fff05288dd3b1b7c32d8..e1cda1fb49681f2f3b4c686bd4f99eb085126032 100644
--- a/src/kalman_likelihood.jl
+++ b/src/kalman_likelihood.jl
@@ -62,7 +62,7 @@ end
KalmanLikelihoodWs(ny, ns, np, nobs) = KalmanLikelihoodWs{Float64, Int64}(ny, ns, np, nobs)
# Z can be either a matrix or a selection vector
-function kalman_likelihood(Y::AbstractArray{U},
+function kalman_likelihood(Y::AbstractArray{Union{U, Missing}},
Z::AbstractArray{W},
H::AbstractArray{U},
T::AbstractArray{U},
@@ -113,7 +113,7 @@ function kalman_likelihood(Y::AbstractArray{U},
return -0.5*(lik_cst + sum(vlik))
end
-function kalman_likelihood(Y::AbstractArray{U},
+function kalman_likelihood(Y::AbstractArray{Union{U, Missing}},
Z::AbstractArray{W},
H::AbstractArray{U},
T::AbstractArray{U},
@@ -244,7 +244,7 @@ function kalman_likelihood_monitored(Y::Matrix{U},
return -0.5*(lik_cst + sum(vlik))
end
-function kalman_likelihood_monitored(Y::AbstractArray{U},
+function kalman_likelihood_monitored(Y::AbstractArray{Union{U, Missing}},
Z::AbstractArray{W},
H::AbstractArray{U},
T::AbstractArray{U},
diff --git a/src/kalman_smoother.jl b/src/kalman_smoother.jl
index 9d260b5f5ae77dc24516a6cefe90a566a086722b..d9df6b36338c98d4a918b284f1702a341e75dd1b 100644
--- a/src/kalman_smoother.jl
+++ b/src/kalman_smoother.jl
@@ -1,4 +1,3 @@
-using Test
struct KalmanSmootherWs{T, U} <: KalmanWs{T, U}
csmall::Vector{T}
Zsmall::Matrix{T}
@@ -30,7 +29,7 @@ struct KalmanSmootherWs{T, U} <: KalmanWs{T, U}
lik::Vector{T}
KT::Matrix{T}
D::Matrix{T}
- ystar::Vector{T}
+ ystar::Vector{Union{T, Missing}}
Zstar::Matrix{T}
Hstar::Matrix{T}
PZi::Vector{T}
@@ -94,7 +93,7 @@ end
KalmanSmootherWs(ny, ns, np, nobs) = KalmanSmootherWs{Float64, Int64}(ny, ns, np, nobs)
-function kalman_smoother!(Y::AbstractArray{U},
+function kalman_smoother!(Y::AbstractArray{Union{U, Missing}},
c::AbstractArray{U},
Z::AbstractArray{W},
H::AbstractArray{U},
@@ -141,9 +140,9 @@ function kalman_smoother!(Y::AbstractArray{U},
pattern = data_pattern[t]
ndata = length(pattern)
vc = changeC ? view(c, :, t) : view(c, :)
- ws.csmall .= view(vc, pattern)
+ # ws.csmall .= view(vc, pattern)
vZsmall = get_vZsmall(ws.Zsmall, ws.iZsmall, Z, pattern, ndata, ny, t)
- vH = changeH ? view(H, :, :, t) : view(H, :, :)
+ vH = changeH ? view(H, pattern, pattern, t) : view(H, pattern, pattern)
vT = changeT ? view(T, :, :, t) : view(T, :, :)
va = changeA ? view(a, :, t) : view(a, :)
vatt = changeAtt ? view(att, :, t) : view(att, :)
@@ -172,16 +171,20 @@ function kalman_smoother!(Y::AbstractArray{U},
update_N!(ws.N, vZsmall, viFZ, ws.L, ws.N_1, ws.PTmp)
end
if length(epsilonh) > 0
- vepsilonh = view(epsilonh, :, t)
+ vepsilonh = view(epsilonh, pattern, t)
# epsilon_t = H*(iF_t*v_t - KDK_t'*r_t) (DK 4.69)
- get_epsilonh!(vepsilonh, vH, viFv, vKDK, ws.r_1, ws.tmp_ny, ws.tmp_ns)
+ vtmp1 = view(ws.tmp_ny, 1:ndata)
+ get_epsilonh!(vepsilonh, vH, viFv, vKDK, ws.r_1, vtmp1, ws.tmp_ns)
if length(Vepsilon) > 0
- vVepsilon = view(Vepsilon,:,:,t)
+ vVepsilon = view(Vepsilon, pattern, pattern, t)
get_iF!(viF, vcholF)
# D_t = inv(F_t) + KDK_t'*N_t*KDK_t (DK 4.69)
- get_D!(ws.D, viF, vKDK, ws.N_1, ws.tmp_ny_ns)
+ vD = view(ws.D, 1:ndata, 1:ndata)
+ vtmp1 = view(ws.tmp_ny_ns, 1:ndata, :)
+ get_D!(vD, viF, vKDK, ws.N_1, vtmp1)
# Vepsilon_t = H - H*D_t*H (DK 4.69)
- get_Vepsilon!(vVepsilon, vH, ws.D, ws.tmp_ny_ny)
+ vtmp1 = view(ws.tmp_ny_ny, 1:ndata, 1:ndata)
+ get_Vepsilon!(vVepsilon, vH, vD, vtmp1)
end
end
if length(etah) > 0
@@ -259,7 +262,7 @@ struct DiffuseKalmanSmootherWs{T, U} <: KalmanWs{T, U}
D::Matrix{T}
uKinf::Vector{T}
uKstar::Vector{T}
- ystar::Vector{T}
+ ystar::Vector{Union{T, Missing}}
Kinf_Finf::Vector{T}
Zstar::Matrix{T}
Hstar::Matrix{T}
@@ -338,7 +341,7 @@ end
DiffuseKalmanSmootherWs(ny, ns, np, nobs) = DiffuseKalmanSmootherWs{Float64, Int64}(ny, ns, np, nobs)
-function diffuse_kalman_smoother_coda!(Y::AbstractArray{U},
+function diffuse_kalman_smoother_coda!(Y::AbstractArray{Union{U, Missing}},
c::AbstractArray{U},
Z::AbstractArray{W},
H::AbstractArray{U},
@@ -399,9 +402,9 @@ function diffuse_kalman_smoother_coda!(Y::AbstractArray{U},
pattern = data_pattern[t]
ndata = length(pattern)
vc = changeC ? view(c, :, t) : view(c, :)
- ws.csmall .= view(vc, pattern)
+ # ws.csmall .= view(vc, pattern)
vZsmall = get_vZsmall(ws.Zsmall, ws.iZsmall, Z, pattern, ndata, ny, t)
- vH = changeH ? view(H, :, :, t) : view(H, :, :)
+ vH = changeH ? view(H, pattern, pattern, t) : view(H, pattern, pattern)
vT = changeT ? view(T, :, :, t) : view(T, :, :)
va = changeA ? view(a, :, t) : view(a, :)
vatt = changeAtt ? view(att, :, t) : view(att, :)
@@ -441,17 +444,19 @@ function diffuse_kalman_smoother_coda!(Y::AbstractArray{U},
ws)
end
if length(epsilonh) > 0
- vepsilonh = view(epsilonh, :, t)
+ vepsilonh = view(epsilonh, pattern, t)
# epsilon_t = -H_t*KDKinf*r0_t (DK p. 135)
- get_epsilonh!(vepsilonh, vH, vKDKinf, r0_1, ws.tmp_ny)
+ vtmp1 = view(ws.tmp_ny, 1:ndata)
+ get_epsilonh!(vepsilonh, vH, vKDKinf, r0_1, vtmp1)
if length(Vepsilon) > 0
- vVepsilon = view(Vepsilon,:,:,t)
+ vVepsilon = view(Vepsilon, pattern, pattern,t)
vTmp = view(ws.tmp_ny_ns, 1:ndata, :)
# D_t = KDKinf_t'*N0_t*KDKinf_t (DK p. 135)
- get_D!(ws.D, vKDK, N0_1, vTmp)
+ vD = view(ws.D, 1:ndata, 1:ndata)
+ get_D!(vD, vKDK, N0_1, vTmp)
# Vepsilon_t = H - H*D_t*H (DK p. 135)
vTmp = view(ws.tmp_ny_ny, 1:ndata, 1:ndata)
- get_Vepsilon!(vVepsilon, vH, ws.D, ws.tmp_ny_ny)
+ get_Vepsilon!(vVepsilon, vH, vD, vTmp)
end
end
if length(etah) > 0
@@ -567,7 +572,7 @@ function diffuse_kalman_smoother_coda!(Y::AbstractArray{U},
end
end
-function diffuse_kalman_smoother!(Y::AbstractArray{U},
+function diffuse_kalman_smoother!(Y::AbstractArray{Union{U, Missing}},
c::AbstractArray{U},
Z::AbstractArray{W},
H::AbstractArray{U},
@@ -615,7 +620,7 @@ function diffuse_kalman_smoother!(Y::AbstractArray{U},
return -0.5*sum(vlik)
end
-function diffuse_kalman_smoother!(Y::AbstractArray{U},
+function diffuse_kalman_smoother!(Y::AbstractArray{Union{U, Missing}},
c::AbstractArray{U},
Z::AbstractArray{W},
H::AbstractArray{U},
diff --git a/src/univariate_step.jl b/src/univariate_step.jl
index 12204f215938799840315fb3ac328a0afbc6ec2e..0218f08ea61750b826b611c3ba4dba129fe6d885 100644
--- a/src/univariate_step.jl
+++ b/src/univariate_step.jl
@@ -313,21 +313,23 @@ function univariate_step(att, a1, Pinftt, Pinf1, Pstartt, Pstar1, Y, t, c, Z, H,
copy!(Pstartt, Pstar)
ny = size(Y,1)
detLTcholH = 1.0
+ ndata = length(pattern)
+ vystar = view(ws.ystar, 1:ndata)
+ vZstar = view(ws.Zstar, 1:ndata, :)
if !isdiag(H)
- detLTcholH = transformed_measurement!(ws.ystar, ws.Zstar, view(Y, :, t), Z, ws.cholH)
+ detLTcholH = transformed_measurement!(ws.ystar, vZstar, view(Y, :, t), Z, ws.cholH)
H = I(ny)
else
- copy!(ws.ystar, view(Y, :, t))
- copy!(ws.Zstar, Z)
+ copy!(vystar, view(Y, pattern, t))
+ copy!(vZstar, Z)
end
llik = 0.0
l2pi = log(2*pi)
- ndata = length(pattern)
for i=1:ndata
uKinf = view(ws.Kinf, i, :, t)
uKstar = view(ws.K, i, :, t)
- Zi = view(ws.Zstar, pattern[i], :)
- ws.v[i] = get_v!(ws.ystar, c, ws.Zstar, a, pattern[i])
+ Zi = view(vZstar, i, :)
+ ws.v[i] = get_v!(vystar, c, vZstar, a, i)
Fstar = get_Fstar!(Zi, Pstartt, H[i], uKstar)
Finf = get_Finf!(Zi, Pinftt, uKinf)
# Conduct check of rank