Commit 6c4488ce by Michel Juillard

### fixed computation of log likelihood and kalman_filter

parent f02dc2fc
 ... ... @@ -25,7 +25,8 @@ using LinearAlgebra.BLAS export KalmanLikelihoodWs, FastKalmanLikelihoodWs, DiffuseKalmanLikelihoodWs export DiffuseKalmanFilterWs export KalmanSmootherWs, kalman_likelihood, kalman_likelihood_monitored export fast_kalman_likelihood, diffuse_kalman_likelihood, kalman_filter! export fast_kalman_likelihood, diffuse_kalman_likelihood export kalman_filter!, diffuse_kalman_filter! export kalman_smoother! abstract type KalmanWs{T, U} end ... ...
 ... ... @@ -239,6 +239,11 @@ function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{T}, z::AbstractVector{ 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} 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} v .= view(y, pattern, t) ... ... @@ -251,12 +256,21 @@ function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{T}, z::AbstractVector{ end # v = y - c - Z*a -- basic function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{T}, c::AbstractArray{T}, z::AbstractArray{T}, a::AbstractArray{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer} 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} 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} v = Y[i] - c[i] @inbounds @simd for j = 1:length(a) v -= Z[i, j]*a[j] end return v 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} copyto!(v, 1, y, iy, ny) ... ... @@ -264,6 +278,11 @@ function get_v!(v::AbstractArray{T}, y::AbstractVecOrMat{T}, c::AbstractArray{T} 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} 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} v .= view(y, pattern, t) .- view(c, pattern) ... ... @@ -298,6 +317,26 @@ function get_Veta!(Veta, Q, R, N, RQ, tmp) end function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, Z::AbstractArray{T}, pattern::AbstractVector{U}, n::U, ny::U) where {T <: AbstractFloat, U <: Integer} vZsmall = view(Zsmall, 1:n, :) if n == ny copyto!(vZsmall, Z) else vZsmall .= view(Z, pattern, :) end return vZsmall end function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, z::AbstractArray{U}, pattern::AbstractVector{U}, n::U, ny::U) where {T <: AbstractFloat, U <: Integer} vZsmall = view(iZsmall, 1:n) if n == ny copyto!(vZsmall, z) else vZsmall .= view(z, pattern) end return vZsmall end function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, Z::AbstractArray{T}, pattern::AbstractVector{U}, n::U, ny::U, t::U) where {T <: AbstractFloat, U <: Integer} changeZ = ndims(Z) > 2 vZ = changeZ ? view(Z, :, :, t) : view(Z, :, :) vZsmall = view(Zsmall, 1:n, :) ... ... @@ -309,7 +348,7 @@ function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, Z::A return vZsmall end function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, z::AbstractArray{U}, pattern::AbstractVector{U}, n::U, ny::U) where {T <: AbstractFloat, U <: Integer} function get_vZsmall(Zsmall::AbstractMatrix{T}, iZsmall::AbstractVector{U}, z::AbstractArray{U}, pattern::AbstractVector{U}, n::U, ny::U, t::U) where {T <: AbstractFloat, U <: Integer} changeZ = ndims(z) > 1 vz = changeZ ? view(z, :, t) : view(z, :,) vZsmall = view(iZsmall, 1:n) ... ...
 ... ... @@ -33,13 +33,12 @@ function kalman_filter!(Y::AbstractArray{U}, vQ = view(Q, :, :, 1) get_QQ!(ws.QQ, vR, vQ, ws.RQ) l2pi = log(2*pi) fill!(ws.lik, 0.0) t = start steady = false vP = view(P, :, :, 1) copy!(ws.oldP, vP) cholHset = false @inbounds while t <= last while t <= last pattern = data_pattern[t] ndata = length(pattern) ... ... @@ -79,7 +78,7 @@ function kalman_filter!(Y::AbstractArray{U}, get_cholF!(ws.cholH, vvH) cholHset = true end ws.lik[t] = univariate_step!(Y, t, ws.Zsmall, vvH, T, ws.QQ, va, vP, ws.kalman_tol, ws) ws.lik[t] = ndata*l2pi + univariate_step!(Y, t, c, ws.Zsmall, vvH, d, T, ws.QQ, va, vP, ws.kalman_tol, ws, pattern) t += 1 continue end ... ... @@ -108,12 +107,8 @@ function kalman_filter!(Y::AbstractArray{U}, end t += 1 end @inbounds if presample > 0 LIK = -0.5*(sum(view(ws.lik,(presample+1):nobs))) else LIK = -0.5*(sum(ws.lik)) end LIK vlik = view(ws.lik, start + presample:last) return -0.5*sum(vlik) end struct DiffuseKalmanFilterWs{T, U} <: KalmanWs{T, U} ... ... @@ -135,8 +130,9 @@ struct DiffuseKalmanFilterWs{T, U} <: KalmanWs{T, U} ZPstar::Matrix{T} Kinf::Matrix{T} iFZ::Matrix{T} Kstar::Matrix{T} K::Matrix{T} PTmp::Matrix{T} oldP::Matrix{T} uKinf::Vector{T} uKstar::Vector{T} Kinf_Finf::Vector{T} ... ... @@ -166,8 +162,9 @@ struct DiffuseKalmanFilterWs{T, U} <: KalmanWs{T, U} ZPstar = Matrix{T}(undef, ny, ns) Kinf = Matrix{T}(undef, ny, ns) iFZ = Matrix{T}(undef, ny, ns) Kstar = Matrix{T}(undef, ny, ns) K = Matrix{T}(undef, ny, ns) PTmp = Matrix{T}(undef, ns, ns) oldP = Matrix{T}(undef, ns, ns) uKinf = Vector{T}(undef, ns) uKstar = Vector{T}(undef, ns) Kinf_Finf = Vector{T}(undef, ns) ... ... @@ -178,7 +175,7 @@ struct DiffuseKalmanFilterWs{T, U} <: KalmanWs{T, U} lik = zeros(T, nobs) kalman_tol = 1e-12 new(csmall, Zsmall, iZsmall, QQ, RQ, c, v, F, iF, iFv, a1, cholF, cholH, ZP, Fstar, ZPstar, Kinf, iFZ, Kstar, PTmp, uKinf, uKstar, Kinf_Finf, ystar, Zstar, Hstar, ZPstar, Kinf, iFZ, K, PTmp, oldP, uKinf, uKstar, Kinf_Finf, ystar, Zstar, Hstar, PZi, lik, kalman_tol) end end ... ... @@ -220,6 +217,7 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U}, vR = view(R, :, :, 1) vQ = view(Q, :, :, 1) get_QQ!(ws.QQ, vR, vQ, ws.RQ) l2pi = log(2*pi) diffuse_kalman_tol = 1e-8 kalman_tol = 1e-8 cholHset = false ... ... @@ -233,6 +231,9 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U}, vT = changeT ? view(T, :, :, t) : view(T, :, :) vR = changeR ? view(R, :, :, t) : view(R, :, :) vQ = changeQ ? view(Q, :, :, t) : view(Q, :, :) if changeR || changeQ get_QQ!(ws.QQ, vR, vQ, ws.RQ) end va = changeA ? view(a, :, t) : view(a, :) va1 = changeA ? view(a, :, t + 1) : view(a, :) vd = changeD ? view(d, :, t) : view(d, :) ... ... @@ -240,9 +241,6 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U}, vPinf1 = changePinf ? view(Pinf, :, :, t + 1) : view(Pinf, :, :) vPstar = changePstar ? view(Pstar, :, :, t) : view(Pstar, :, :) vPstar1 = changePstar ? view(Pstar, :, :, t + 1) : view(Pstar, :, :) if changeR || changeQ get_QQ!(ws.QQ, vR, vQ, ws.RQ) end vv = view(ws.v, 1:ndata) vvH = view(vH, pattern, pattern) vZP = view(ws.ZP, 1:ndata, :) ... ... @@ -256,7 +254,7 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U}, vcholH = view(ws.cholH, 1:ndata, 1:ndata, 1) viFv = view(ws.iFv, 1:ndata) vKinf = view(ws.Kinf, 1:ndata, :, 1) vKstar = view(ws.Kstar, 1:ndata, :, 1) vKstar = view(ws.K, 1:ndata, :, 1) # v = Y[:,t] - c - Z*a get_v!(vv, Y, vc, vZsmall, va, t, pattern) ... ... @@ -271,10 +269,10 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U}, get_cholF!(vcholH, vvH) cholHset = true end ws.lik[t] += univariate_step(Y, t, vZsmall, vvH, T, QQ, a, Pinf, Pstar, diffuse_kalman_tol, kalman_tol, ws) ws.lik[t] += ndata*l2pi + univariate_step(Y, t, c, vZsmall, vvH, d, T, QQ, a, Pinf, Pstar, diffuse_kalman_tol, kalman_tol, ws, pattern) end else ws.lik[t] = log(det_from_cholesky(vcholF)) ws.lik[t] = ndata*l2pi + log(det_from_cholesky(vcholF)) # Kinf = iFinf*Z*Pinf %define Kinf'=T^{-1}*K_0 with M_{\infty}=Pinf*Z' copy!(vKinf, vZP) LAPACK.potrs!('U', vcholF, vKinf) ... ... @@ -296,66 +294,55 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U}, t end function diffuse_kalman_filter(Y::Matrix{U}, Z::AbstractArray{W}, H::Matrix{U}, T::Matrix{U}, R::Matrix{U}, Q::Matrix{U}, a::Vector{U}, Pinf::Matrix{U}, Pstar::Matrix{U}, start::V, last::V, presample::V, tol::U, ws::DiffuseKalmanLikelihoodWs) where {U <: AbstractFloat, V <: Integer, W <: Real} function diffuse_kalman_filter!(Y::AbstractArray{U}, c::AbstractArray{U}, Z::AbstractArray{W}, H::AbstractArray{U}, d::AbstractArray{U}, T::AbstractArray{U}, R::AbstractArray{U}, Q::AbstractArray{U}, a::AbstractArray{U}, Pinf::AbstractArray{U}, Pstar::AbstractArray{U}, start::V, last::V, presample::V, tol::U, ws::DiffuseKalmanFilterWs, data_pattern::Vector{Vector{V}}) where {U <: AbstractFloat, V <: Integer, W <: Real} ny = size(Y,1) nobs = last - start + 1 get_QQ!(ws.QQ, R, Q, ws.RQ) lik_cst = (nobs - presample)*ny*log(2*pi) fill!(ws.lik, 0.0) t = diffuse_kalman_likelihood_init!(Y, Z, H, T, ws.QQ, a, Pinf, Pstar, start, last, tol, ws) kalman_likelihood(Y, Z, H, T, R, Q, a, Pstar, t, last, presample, ws) @inbounds if presample > 0 LIK = -0.5*(lik_cst + sum(view(ws.lik, (presample+1):nobs))) else LIK = -0.5*(lik_cst + sum(ws.lik)) end LIK t = diffuse_kalman_filter_init!(Y, c, Z, H, d, T, R, Q, a, Pinf, Pstar, start, last, presample, tol, ws, data_pattern) kalman_filter!(Y, c, Z, H, d, T, R, Q, a, Pstar, t + 1, last, presample, ws, data_pattern) vlik = view(ws.lik, start + presample:last) return -0.5*sum(vlik) end function diffuse_kalman_filter!(Y::Matrix{U}, Z::AbstractArray{W}, H::Matrix{U}, T::Matrix{U}, R::Matrix{U}, Q::Matrix{U}, a::Vector{U}, Pinf::Matrix{U}, Pstar::Matrix{U}, start::V, last::V, presample::V, tol::U, ws::DiffuseKalmanFilterWs, data_pattern::Vector{Vector{V}}) where {U <: AbstractFloat, V <: Integer, W <: Real} ny = size(Y,1) nobs = last - start + 1 get_QQ!(ws.QQ, R, Q, ws.RQ) lik_cst = (nobs - presample)*ny*log(2*pi) fill!(ws.lik, 0.0) t = diffuse_kalman_filter_init!(Y, Z, H, T, ws.QQ, a, Pinf, Pstar, start, last, tol, ws, data_pattern) kalman_likelihood(Y, Z, H, T, R, Q, a, Pstar, t, last, presample, ws, data_pattern) @inbounds if presample > 0 LIK = -0.5*(lik_cst + sum(view(ws.lik, (presample+1):nobs))) else LIK = -0.5*(lik_cst + sum(ws.lik)) end LIK end function diffuse_kalman_filter!(Y::AbstractArray{U}, c::AbstractArray{U}, Z::AbstractArray{W}, H::AbstractArray{U}, d::AbstractArray{U}, T::AbstractArray{U}, R::AbstractArray{U}, Q::AbstractArray{U}, a::AbstractArray{U}, Pinf::AbstractArray{U}, Pstar::AbstractArray{U}, start::V, last::V, presample::V, tol::U, ws::DiffuseKalmanFilterWs) where {U <: AbstractFloat, V <: Integer, W <: Real} m, n = size(Y) full_data_pattern = [collect(1:m) for i = 1:n] diffuse_kalman_filter!(Y, c, Z, H, d, T, R, Q, a, Pinf, Pstar, start, last, presample, tol, ws, full_data_pattern) end
 ... ... @@ -78,8 +78,6 @@ function kalman_likelihood(Y::AbstractArray{U}, nobs = last - start + 1 # QQ = R*Q*R' get_QQ!(ws.QQ, R, Q, ws.RQ) lik_cst = (nobs - presample)*ny*log(2*pi) fill!(ws.lik, 0.0) cholHset = false t = start iy = (start - 1)*ny + 1 ... ... @@ -110,12 +108,9 @@ function kalman_likelihood(Y::AbstractArray{U}, end t += 1 end @inbounds if presample > 0 LIK = -0.5*(lik_cst + sum(view(ws.lik,(presample+1):nobs))) else LIK = -0.5*(lik_cst + sum(ws.lik)) end LIK lik_cst = (nobs - presample)*ny*log(2*pi) vlik = view(ws.lik, start + presample:last) return -0.5*(lik_cst + sum(vlik)) end function kalman_likelihood(Y::AbstractArray{U}, ... ... @@ -136,7 +131,6 @@ function kalman_likelihood(Y::AbstractArray{U}, # QQ = R*Q*R' get_QQ!(ws.QQ, R, Q, ws.RQ) l2pi = log(2*pi) fill!(ws.lik, 0.0) t = start iy = (start - 1)*ny + 1 ncolZ = size(Z, 2) ... ... @@ -161,13 +155,12 @@ function kalman_likelihood(Y::AbstractArray{U}, get_F!(vF, vZP, vZsmall, P, vH) info = get_cholF!(vcholF, vF) if info != 0 @show info # F is near singular if !cholHset get_cholF!(vcholH, vH) cholHset = true end ws.lik[t] = univariate_step!(Y, t, Z, H, T, ws.QQ, a, P, ws.kalman_tol, ws) ws.lik[t] = ndata*l2pi + univariate_step!(Y, t, Z, H, T, ws.QQ, a, P, ws.kalman_tol, ws, pattern) else # iFv = inv(F)*v get_iFv!(viFv, vcholF, vv) ... ... @@ -181,12 +174,8 @@ function kalman_likelihood(Y::AbstractArray{U}, end t += 1 end @inbounds if presample > 0 LIK = -0.5*sum(view(ws.lik,(presample+1):nobs)) else LIK = -0.5*sum(ws.lik) end LIK vlik = view(ws.lik, start + presample:last) return -0.5*sum(vlik) end function kalman_likelihood_monitored(Y::Matrix{U}, ... ... @@ -206,8 +195,6 @@ function kalman_likelihood_monitored(Y::Matrix{U}, ns = size(T,1) # QQ = R*Q*R' get_QQ!(ws.QQ, R, Q, ws.RQ) lik_cst = (nobs - presample)*ny*log(2*pi) fill!(ws.lik, 0.0) t = start iy = (start - 1)*ny + 1 steady = false ... ... @@ -221,7 +208,6 @@ function kalman_likelihood_monitored(Y::Matrix{U}, get_F!(ws.F, ws.ZP, Z, P, H) info = get_cholF!(ws.cholF, ws.F) if info != 0 @show info # F is near singular if !cholHset get_cholF!(ws.cholH, H) ... ... @@ -253,12 +239,9 @@ function kalman_likelihood_monitored(Y::Matrix{U}, end t += 1 end @inbounds if presample > 0 LIK = -0.5*(lik_cst + sum(view(ws.lik,(presample+1):nobs))) else LIK = -0.5*(lik_cst + sum(ws.lik)) end LIK lik_cst = (nobs - presample)*ny*log(2*pi) vlik = view(ws.lik, start + presample:last) return -0.5*(lik_cst + sum(vlik)) end function kalman_likelihood_monitored(Y::AbstractArray{U}, ... ... @@ -280,7 +263,6 @@ function kalman_likelihood_monitored(Y::AbstractArray{U}, # QQ = R*Q*R' get_QQ!(ws.QQ, R, Q, ws.RQ) l2pi = log(2*pi) fill!(ws.lik, 0.0) t = start iy = (start - 1)*ny + 1 steady = false ... ... @@ -310,7 +292,7 @@ function kalman_likelihood_monitored(Y::AbstractArray{U}, get_cholF!(ws.cholH, H) cholHset = true end ws.lik[t] = univariate_step!(Y, t, Z, H, T, ws.QQ, a, P, ws.kalman_tol, ws) ws.lik[t] = ndatat*l2pi + univariate_step!(Y, t, Z, H, T, ws.QQ, a, P, ws.kalman_tol, ws, pattern) t += 1 continue end ... ... @@ -336,12 +318,8 @@ function kalman_likelihood_monitored(Y::AbstractArray{U}, end t += 1 end @inbounds if presample > 0 LIK = -0.5*(sum(view(ws.lik,(presample+1):nobs))) else LIK = -0.5*(sum(ws.lik)) end LIK vlik = view(ws.lik, start + presample:last) return -0.5*sum(vlik) end struct FastKalmanLikelihoodWs{T, U} <: KalmanWs{T, U} ... ... @@ -425,8 +403,6 @@ function fast_kalman_likelihood(Y::Matrix{U}, ny = size(Y,1) nobs = last - start + 1 get_QQ!(ws.QQ, R, Q, ws.RQ) lik_cst = (nobs - presample)*ny*log(2*pi) fill!(ws.lik, 0.0) # F = Z*P*Z' + H get_F!(ws.F, ws.ZP, Z, P, H) get_cholF!(ws.cholF, ws.F) ... ... @@ -460,12 +436,9 @@ function fast_kalman_likelihood(Y::Matrix{U}, update_W!(ws.W, ws.ZW, ws.cholF, T, ws.K, ws.iFZW, ws.KtiFZW) t += 1 end @inbounds if presample > 0 LIK = -0.5*(lik_cst + sum(view(ws.lik, (presample+1):nobs))) else LIK = -0.5*(lik_cst + sum(ws.lik)) end LIK lik_cst = (nobs - presample)*ny*log(2*pi) vlik = view(ws.lik, start + presample:last) return -0.5*(lik_cst + sum(vlik)) end function fast_kalman_likelihood(Y::Matrix{U}, ... ... @@ -485,7 +458,6 @@ function fast_kalman_likelihood(Y::Matrix{U}, nobs = last - start + 1 get_QQ!(ws.QQ, R, Q, ws.RQ) l2pi = log(2*pi) fill!(ws.lik, 0.0) # F = Z*P*Z' + H get_F!(ws.F, ws.ZP, Z, P, H) get_cholF!(ws.cholF, ws.F) ... ... @@ -534,12 +506,8 @@ function fast_kalman_likelihood(Y::Matrix{U}, update_W!(vW, vZW, vcholF, T, vK, viFZW, vKtiFZW) t += 1 end @inbounds if presample > 0 LIK = -0.5*(sum(view(ws.lik, (presample+1):nobs))) else LIK = -0.5*(sum(ws.lik)) end LIK vlik = view(ws.lik, start + presample:last) return -0.5*sum(vlik) end struct DiffuseKalmanLikelihoodWs{T, U} <: KalmanWs{T, U} ... ... @@ -686,7 +654,6 @@ function diffuse_kalman_likelihood_init!(Y::Matrix{U}, ny = size(Y, 1) t = start LIK = 0 iy = (start - 1)*ny + 1 diffuse_kalman_tol = 1e-8 kalman_tol = 1e-8 ... ... @@ -715,7 +682,7 @@ function diffuse_kalman_likelihood_init!(Y::Matrix{U}, if norm(vF) < tol return t - 1 else ws.lik[t] += univariate_step(Y, t, vZsmall, H, T, QQ, a, Pinf, Pstar, diffuse_kalman_tol, kalman_tol, pattern, ws) ws.lik[t] += ndata*l2pi + univariate_step(Y, t, vZsmall, H, T, QQ, a, Pinf, Pstar, diffuse_kalman_tol, kalman_tol, pattern, ws, pattern) end else ws.lik[t] = log(det_from_cholesky(ws.cholF)) ... ... @@ -765,16 +732,11 @@ function diffuse_kalman_likelihood(Y::Matrix{U}, ny = size(Y,1) nobs = last - start + 1 get_QQ!(ws.QQ, R, Q, ws.RQ) lik_cst = (nobs - presample)*ny*log(2*pi) fill!(ws.lik, 0.0) t = diffuse_kalman_likelihood_init!(Y, Z, H, T, ws.QQ, a, Pinf, Pstar, start, last, tol, ws) kalman_likelihood(Y, Z, H, T, R, Q, a, Pstar, t, last, presample, ws) @inbounds if presample > 0 LIK = -0.5*(lik_cst + sum(view(ws.lik, (presample+1):nobs))) else LIK = -0.5*(lik_cst + sum(ws.lik)) end LIK kalman_likelihood(Y, Z, H, T, R, Q, a, Pstar, t + 1, last, presample, ws) lik_cst = (nobs - presample)*ny*log(2*pi) vlik = view(ws.lik, start + presample:last) return -0.5*(lik_cst + sum(vlik)) end function diffuse_kalman_likelihood(Y::Matrix{U}, ... ... @@ -797,15 +759,9 @@ function diffuse_kalman_likelihood(Y::Matrix{U}, ny = size(Y,1) nobs = last - start + 1 get_QQ!(ws.QQ, R, Q, ws.RQ) lik_cst = (nobs - presample)*ny*log(2*pi) fill!(ws.lik, 0.0) t = diffuse_kalman_likelihood_init!(Y, Z, H, T, ws.QQ, a, Pinf, Pstar, start, last, tol, ws, data_pattern) kalman_likelihood(Y, Z, H, T, R, Q, a, Pstar, t, last, presample, ws, data_pattern) @inbounds if presample > 0 LIK = -0.5*(lik_cst + sum(view(ws.lik, (presample+1):nobs))) else LIK = -0.5*(lik_cst + sum(ws.lik)) end LIK kalman_likelihood(Y, Z, H, T, R, Q, a, Pstar, t + 1, last, presample, ws, data_pattern) vlik = view(ws.lik, start + presample:last) return -0.5*sum(vlik) end
 ... ... @@ -54,6 +54,38 @@ function univariate_step!(Y, t, Z, H, T, RQR, a, P, kalman_tol, ws) return llik + 2*log(detLTcholH) end function univariate_step!(Y, t, Z, H, T, RQR, a, P, kalman_tol, ws, pattern) ny = size(Y,1) detLTcholH = 1.0 if !isdiag(H) detLTcholH = transformed_measurement!(ws.ystar, ws.Zstar, view(Y, :, t), Z, ws.cholH) H = I(ny) else copy!(ws.ystar, view(Y, :, t)) copy!(ws.Zstar, Z) end l2pi = log(2*pi) llik = 0.0 ndata = size(pattern) for i=1:ndata Zi = view(ws.Zstar, pattern[i], :) v = get_v!(ws.ystar, ws.Zstar, a, pattern[i]) F = get_F(Zi, P, H[i,i], ws.PZi) if abs(F) > kalman_tol a .+= (v/F) .* ws.PZi # P = P - PZi*PZi'/F ger!(-1.0/F, ws.PZi, ws.PZi, P) llik += ndata*l2pi + log(F) + v*v/F end end mul!(ws.a1, T, a) a .= ws.a1 mul!(ws.PTmp, T, P) copy!(P, RQR) mul!(P, ws.PTmp, T', 1.0, 1.0) return llik + 2*log(detLTcholH) end function univariate_step(Y, t, Z, H, T, QQ, a, Pinf, Pstar, diffuse_kalman_tol, kalman_tol, ws) ny = size(Y,1) detLTcholH = 1.0 ... ... @@ -67,6 +99,7 @@ function univariate_step(Y, t, Z, H, T, QQ, a, Pinf, Pstar, diffuse_kalman_tol, llik = 0.0