%% Kalman filter

n = 100; %number of trials to simulate
sh = 1; %variance of random walk step
se = 10; %variance of observation noise
m = zeros(n+1,1); %sequence of prior means; zero for first trial
s = zeros(n+1,1); s(1) = sh; %sequence of prior variances; sh for first trial

x = cumsum(randn(n,1)*sqrt(sh)); %sequence of true values
y = x + randn(n,1)*sqrt(se); %sequence of observations

for i=1:n
    m(i+1) = (m(i)*se+y(i)*s(i))/(se+s(i));
    s(i+1) = se*s(i)/(se+s(i)) + sh;
end

figure(4),clf,hold on
plot(x,'k')
plot(y,'r')
plot(m,'b')
b = fill([1:n+1,n+1:-1:1],[m-sqrt(s);flipud(m+sqrt(s))],[.9 .9 .9],'linestyle','none');
alpha(b,.5)

%% Markov chain stationary distribution

T = [.6 .1 .1;.3 .8 0;.1 .1 .9]; %transition matrix
p = [.2;.3;.5]; %stationary distribution

disp('Trasition matrix (T):'),disp(T)
disp('Stationary distributions (p):'),disp(p)
disp('T*p:'),disp(T*p)

n = 10000; %number of steps to run
burn = 1000; %initial steps to discard
space = 10; %spacing of samples to save

sample = zeros((n-burn)/space,1);
x = randi(3); %random starting state
for i=1:n
    x = find(rand<cumsum(T(:,x)),1); %sample according to xth column of T
    if i>burn && mod(i,space)==0
        sample((i-burn)/space) = x; %save this sample
    end
end

p = [mean(sample==1) mean(sample==2) mean(sample==3)]; %proportion for each state
disp(['Simulated proportions: ' num2str(p)]) %display proportions
disp(['Standard errors: ' num2str(sqrt(p.*(1-p)*space/(n-burn)))]) %standard errors assuming independence
