%Constructing data-driven Markov Chains with Matlab
%Example 3.
%(JESSE DORRESTIJN, 7 June 2016)
% Let us now introduce conditioning. Imagine that the transition
% probabilities depend on a certain variable X. Then, the Markov chain
% becomes a Conditional Markov Chain, because it is conditioned on X. The
% observational data now consists of a sequence y_obs (as was the case in
% example 1 and 2) and an additional sequence X_obs. After construction of
% the Conditional Markov Chain, an additional sequence X should be
% available, e.g. X=X_obs.
%
clear
clc
close all
%YOUR INPUT:
M = 5; %Number of states of the sequence X.
N = 5; %Number of states of the Markov chain.
L = 100000; %length of observational sequence.
T = 1000; %Plot the first T states of the sequence. T should be smaller than L.
%First we produce two sequence that we later use as observational data
%for the construction of the Markov chain: y_obs and X_obs.
y_obs = zeros(L,1); %y_obs will be the input or OBServational sequence.
PX = zeros(N,N,M);
PX_cum = PX;
for k = 1:M %For each of the M states of X, we produce a sequence with
P = rand(N); %transition matrix used to construct input sequence.
for j=1:N
P(:,j) = P(:,j)*(N-j+1)^((M-k)*(M-k));
end
P_cum = P;
for j=2:N
P_cum(:,j) = P_cum(:,j-1) + P(:,j); %cumulative version of P.
end
for j=1:N
P(:,j) = P(:,j)./P_cum(:,N); %Normalize transition matrix
end
P_cum = P;
for j=2:N
P_cum(:,j) = P_cum(:,j-1) + P(:,j); %normalized cumulative version of P.
end
PX(:,:,k) = P;
PX_cum(:,:,k) = P_cum;
end
%Construct sequence for X:
Q = rand(M); %transition matrix used to construct input sequence.
Q = Q+Q.*eye(M)*100;
Q_cum = Q;
for j=2:M
Q_cum(:,j) = Q_cum(:,j-1) + Q(:,j); %cumulative version of P.
end
for j=1:M
Q(:,j) = Q(:,j)./Q_cum(:,M); %Normalize transition matrix
end
Q_cum = Q;
for j=2:M
Q_cum(:,j) = Q_cum(:,j-1) + Q(:,j); %normalized cumulative version of P.
end
X_obs = zeros(L,1);
X_obs(1) = 1; %sequence starts with a 1.
for t=1:L-1 %construct entire sequence X_obs
r = rand;
X_obs(t+1) = sum(r>Q_cum(X_obs(t),:))+1;
end
y_obs(1) = 1; %sequence starts with a 1;
for t=1:L-1 %construct entire sequence y_obs
r = rand;
y_obs(t+1) = sum(r>PX_cum(y_obs(t),:,X_obs(t)))+1;
end
%Plot the sequence of the observations:
figure(1)
subplot(3,1,1)
plot(y_obs(1:T)) %plot the first T states of the sequence y_obs
hold on
plot(X_obs(1:T),'r--'); %plot the first T states of the sequence X_obs
xlim([-10 T+10])
ylim([0.5 N+0.5])
set(gca,'YTick',1:N)
xlabel('observations')
ylabel('state')
%%
%Estimation of transition probability matrix of Markov Chain.
P_MC = zeros(N,N);
for t=1:L-1
P_MC(y_obs(t),y_obs(t+1))= P_MC(y_obs(t),y_obs(t+1))+1;
end
P_MC_cum = P_MC;
for j=2:N
P_MC_cum(:,j) = P_MC_cum(:,j-1) + P_MC(:,j); %cumulative version of P.
end
for j=1:N
P_MC(:,j) = P_MC(:,j)./P_MC_cum(:,N); %Normalize transition matrix
end
P_MC_cum = P_MC;
for j=2:N
P_MC_cum(:,j) = P_MC_cum(:,j-1) + P_MC(:,j); %normalized cumulative version of P.
end
%%
%Produce sequence with the Markov chain.
y_MC = zeros(L,1); %y_obs will be the input or OBServational sequence.
y_MC(1) = 1; %sequence starts with a 1;
for t=1:L-1 %construct entire sequence
r = rand;
y_MC(t+1) = sum(r>P_MC_cum(y_MC(t),:))+1;
end
%Plot sequence:
figure(1)
subplot(3,1,2)
plot(y_MC(1:T)) %plot the first T states of the sequence (T defined above).
xlim([-10 T+10])
ylim([0.5 N+0.5])
set(gca,'YTick',1:N)
xlabel('Markov Chain')
ylabel('state')
%%
%Estimation of transition probability matrices of Conditional Markov Chain.
PX_CMC = zeros(N,N,M);
for t=1:L-1
PX_CMC(y_obs(t),y_obs(t+1),X_obs(t))= PX_CMC(y_obs(t),y_obs(t+1),X_obs(t))+1;
end
PX_CMC_cum = PX_CMC;
for k=1:M
for j=2:N
PX_CMC_cum(:,j,k) = PX_CMC_cum(:,j-1,k) + PX_CMC(:,j,k); %cumulative version of PX_MC.
end
for j=1:N
PX_CMC(:,j,k) = PX_CMC(:,j,k)./PX_CMC_cum(:,N,k); %Normalize transition matrix
end
PX_CMC_cum(:,:,k) = PX_CMC(:,:,k);
for j=2:N
PX_CMC_cum(:,j,k) = PX_CMC_cum(:,j-1,k) + PX_CMC(:,j,k); %normalized cumulative version of P.
end
end
%%
%Produce sequence with the Conditional Markov chain.
y_CMC = zeros(L,1); %y_obs will be the input or OBServational sequence.
y_CMC(1) = 1; %sequence starts with a 1;
for t=1:L-1 %construct entire sequence
r = rand;
y_CMC(t+1) = sum(r>PX_CMC_cum(y_CMC(t),:,X_obs(t)))+1;
end
%Plot sequence:
figure(1)
subplot(3,1,3)
plot(y_CMC(1:T)) %plot the first T states of the sequence (T defined above).
hold on
plot(X_obs(1:T),'r--')
xlim([-10 T+10])
ylim([0.5 N+0.5])
set(gca,'YTick',1:N)
xlabel('Conditional Markov Chain')
ylabel('state')