From 934f8ef59ae0e771b0af97d3db9f7b7b3c07b5ee Mon Sep 17 00:00:00 2001 From: MichelJuillard Date: Mon, 13 Feb 2023 15:52:42 +0100 Subject: [PATCH] WIP univariate diffuse smoother --- doc/notes.tex | 3 +-- src/univariate_step.jl | 7 +++++++ test/test_univariate_step.jl | 8 +++++--- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/doc/notes.tex b/doc/notes.tex index 321a66b..5c4e083 100644 --- a/doc/notes.tex +++ b/doc/notes.tex @@ -236,8 +236,7 @@ \section{Univariate diffuse smoother step} \end{align*} For $i = p_{t-1},\ldots,0$, if $|F_{t,i}| > 0$, -\begin{align*} - +\begin{align*} r0_{t-1,i-1} &= L'_{\infty,i}r0_{t,i} \\ r1_{t-1,i-1} &= Z'_{t,i}F^{-1}_{t,i}\nu_{t,i} + L'_{\infty,t,i}r0_{t,i} + L'_{0,t,i}r1_{t,i} \\ diff --git a/src/univariate_step.jl b/src/univariate_step.jl index f231254..e2c5579 100644 --- a/src/univariate_step.jl +++ b/src/univariate_step.jl @@ -223,6 +223,7 @@ function diffuse_univariate_step!(Y, t, Z, H, T, RQR, a, Pinf, Pstar, diffuse_ka if Finf > diffuse_kalman_tol # F_{\infty,t,i} = 0, use upper part of bracket on p. 175 DK (2012) for w_{t,i} copy!(K0, ZPinf) rmul!(K0, 1/Finf) + @show K0 a .+= v.*K0 # Pstar = Pstar + K0*(K0_Finf'*(Fstar/Finf)) - K1*K0_Finf' - K0_Finf*K1' ger!( Fstar/Finf, ZPinf, K0, Pstar) @@ -401,6 +402,7 @@ function univariate_diffuse_smoother_step!(T::AbstractMatrix{W}, ws::DiffuseKalmanSmootherWs) where {W <: AbstractFloat, U <: Real} ny = size(Finf, 1) ns = size(L0, 1) + @show v @show r0_1 @show r1_1 for i = ny: -1: 1 @@ -418,12 +420,17 @@ function univariate_diffuse_smoother_step!(T::AbstractMatrix{W}, # get_L1!(L1, K0, K1, Finf, Fstar, Z, i) fill!(L1, 0.0) ger!(-1.0, vK1, vZ, L1) + @show vK1 + @show vZ # compute r1_{t,i-1} first because it depends # upon r0{t,i} # update_r1!(r1, r0, Z, v, Finv, L0, L1, i) copy!(ws.tmp_ns, vZ) rmul!(ws.tmp_ns, v[i]*iFinf) + @show ws.tmp_ns + @show L1 mul!(ws.tmp_ns, transpose(L1), r0_1, 1.0, 1.0) + @show ws.tmp_ns mul!(ws.tmp_ns, transpose(L0), r1_1, 1.0, 1.0) copy!(r1, ws.tmp_ns) # r0(:,t) = Linf'*r0(:,t); % KD (2000), eq. (25) for r_0 diff --git a/test/test_univariate_step.jl b/test/test_univariate_step.jl index 8c55249..a6a142c 100644 --- a/test/test_univariate_step.jl +++ b/test/test_univariate_step.jl @@ -380,6 +380,8 @@ pstar0 = copy(Pstar[:, :, 1]) pstar = copy(pstar0) QQ = R*Q*R' KalmanFilterTools.diffuse_univariate_step!(y, t, Z, H, T, QQ, a, pinf, pstar, 1e-10, 1e-10, ws_d) +@show ws_d.K0[:, :, 1] +@show ws_d.K[:, :, 1] v = y - c - Z*a0 K0 = T*pinf0*Z'*inv(Z*pinf0*Z') a1 = T*a0 + K0*v @@ -396,8 +398,8 @@ K = copy(K0) L0_target = I(ns) - K0*Z L1_target = -K1*Z -r1_target = Z'inv(Finf)*v + L0_target'*r1 +L1_target'*r0 -r0_target = L0_target'*r0 +r1_target = Z'inv(Finf)*v + L0_target'*r1_1 +L1_target'*r0_1 +r0_target = L0_target'*r0_1 N0_target = L0_target'*N0*L0_target N1_target = Z'*inv(Finf)*Z + L0_target'*N1*L0_target + L1_target'*N0_target*L0_target Fstar = -inv(Finf)*Fstar*inv(Finf) @@ -407,7 +409,7 @@ tol = 1e-12 KalmanFilterTools.univariate_diffuse_smoother_step!(T, ws_d.F[:, :, 1], ws_d.Fstar[:, :, 1], ws_d.K0[:, :, 1], ws_d.K[:, :, 1], ws_d.L, ws_d.L1, ws_d.N, ws_d.N1, -ws_d.N2, r0, r1, ws_d.v[:,1], Z, +ws_d.N2, r0, r1, r0_1, r1_1, ws_d.v[:,1], Z, pinf, pstar, tol, ws_d) #=