Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- function [x_T, P_T, P_Tm2] = smoother(G,F,Q,R,mu0,Sigma,y)
- % G is T*p
- p = size(mu0,1); % mu0 is p*1
- T = size(y,1); % y is T*1
- J = zeros(p,p,T);
- K = zeros(p,T); % gain matrix
- x = zeros(p,T);
- x(:,1) = mu0;
- x_m1 = zeros(p,T);
- x_T = zeros(p,T); % x values when we know all the data
- % Notation : x = xt given t ; x_m1 = xt given t-1 (m1 stands for minus
- % one)
- P = zeros(p,p,T);% array of cov(xt|y1...yt), eq (6) in Shumway & Stoffer 1982
- P(:,:,1) = Sigma;
- P_m1 = zeros(p,p,T); % Same notation ; = cov(xt, xt-1|y1...yt) , eq (7)
- P_T = zeros(p,p,T);
- P_Tm2 = zeros(p,p,T); % cov(xT, xT-1|y1...yT)
- for t = 2:T %starts at t = 2 because at each time t we need info about t-1
- x_m1(:,t) = F*x(:,t-1); % eq A3 ; pxp * px1 = px1
- P_m1(:,:,t) = (F*P(:,:,t-1)*F') + Q; % A4 ; pxp * pxp = pxp
- if nnz(isnan(P_m1(:,:,t)))
- error('NaNs in P_m1 at time t = %d',t);
- end
- if nnz(isinf(P_m1(:,:,t)))
- error('Infs in P_m1 at time t = %d',t);
- end
- K(:,t) = P_m1(:,:,t)*G(t,:)'*pinv((G(t,:)*P_m1(:,:,t)*G(t,:)') + R); %A5 ; pxp * px1 * 1*1 = p*1
- %K(:,t) = P_m1(:,:,t)*G(t,:)'/((G(t,:)*P_m1(:,:,t)*G(t,:)') + R); %A5 ; pxp * px1 * 1*1 = p*1
- % The matrix inversion seems to generate NaN values which quickly
- % contaminate all the other matrices. There is no warning about
- % (close to) singular matrices or whatever. The use of pinv()
- % instead of inv() seems to solve the problem... but I don't think
- % it's the appropriate way to deal with it, there must be something
- % wrong elsewhere
- if nnz(isnan(K(:,t)))
- error('NaNs in K at time t = %d',t);
- end
- x(:,t) = x_m1(:,t) + (K(:,t)*(y(t)-(G(t,:)*x_m1(:,t)))); %A6
- P(:,:,t) = P_m1(:,:,t) - (K(:,t)*G(t,:)*P_m1(:,:,t)); %A7
- end
- x_T(:,T) = x(:,T);
- P_T(:,:,T) = P(:,:,T);
- for t = T:-1:2 % we stop at 2 since we need to use t-1
- %P_m1 seem to get really huge (x10^22...), might lead to "Inf"
- %values which in turn might screw pinv()
- %% inv() caused NaN value to appear, pinv seems to solve the issue
- J(:,:,t-1) = P(:,:,t-1)*F'*pinv(P_m1(:,:,t)); % A8 pxp * pxp * pxp
- %J(:,:,t-1) = P(:,:,t-1)*F'/(P_m1(:,:,t)); % A8 pxp * pxp * pxp
- x_T(:,t-1) = x(:,t-1) + J(:,:,t-1)*(x_T(:,t)-(F*x(:,t-1))); %A9 % Becomes NaN during 8th iteration!
- P_T(:,:,t-1) = P(:,:,t-1) + J(:,:,t-1)*(P_T(:,:,t)-P_m1(:,:,t))*J(:,:,t-1)'; %A10
- nans = [nnz(isnan(J)) nnz(isnan(P_m1)) nnz(isnan(F)) nnz(isnan(x_T)) nnz(isnan(x_m1))];
- if nnz(nans)
- error('NaN invasion at time t = %d',t);
- end
- end
- P_Tm2(:,:,T) = (eye(p) - K(:,T)*G(T,:))*F*P(:,:,T-1); % %A12
- for t = T:-1:3 % stop at 3 because use of t-2
- P_Tm2(:,:,t-1) = P_m1(:,:,t-1)*J(:,:,t-2)' + J(:,:,t-1)*(P_Tm2(:,:,t)-F*P(:,:,t-1))*J(:,:,t-2)'; % A11
- end
- end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement