clear; set(0,'defaultaxesfontsize',14); format long
%%% smoothing problem for a deterministic scalar OU process
%%% based on Ch. 2.8 in: Law, Stuart, Zygalakis (2015), "Data Assimilation"


N=300% number of steps
A=-0.2;% dynamics determined by 
m0 = 4; % prior mean
sigma0 = sqrt(5); % prior variance is sigma0^2; try different values between 1 and 0.01
gamma=0.1;% observation noise variance is gamma^25
sd=1;rng(sd);% choose random number seed; try different seeds, e.g. 0,1,2,3,...
x=zeros(N+1,1);  y=zeros(N+1,1);% preallocate space
x(1,:)=1; % initial condition

%% generate data
for i=1:N
    x(i+1,:)=x(i,:)*A';
    y(i+1,:)=x(i+1,:)+gamma.*randn(1);
end

figure(1)
plot(x(:,1)); xlim([0 20])
hold on 
plot(y(2:N+1,1),'o')
hold off
title('dynamics')



%% posterior mean and variance

LAMBDA = A.^(1:N);

sigma1 = (1/sigma0^2 + A^2/gamma^2 * (1-A^(2*N)/(1-A^2)))^(-1/2)
m1 = sigma1^2 * (m0/(sigma0^2) + 1/(gamma^2) * (LAMBDA*y(2:N+1,1)) )

t = [-10:0.1:10];

figure(2)
plot(t,normpdf(t,m0,sigma0))
hold on
plot(t,normpdf(t,m1,sigma1))
xline(x(1,1),':')
hold off
title('smoothing distribution')
legend('prior','posterior','truth')

