function Tuesday_project;


%
% An example of first order Markov chain. 
% p(x_i=l|x_{i-1}=k,x_{i-2}=k2,....)=p(x_i=l|x_{i-1}=k)
%
% 1. Compute marginal probabilities of success ( and draw realizations from markov chains).
%
% Compute the VOI for imperfect information at data_nodes
%
% PV=sum_i max (-CostV'*p_i(x_i),-Cost);
% PoV(x_j)= sum_{x_j} sum_i max (-CostV'*p_i(x_i|x_j),-Cost) * p(x_j)
% VOI(x_j) = PoV(x_j)-PV 
%
% Subfile (at bottom of this file: evidence_prop_Markov_chain)
% - computing p(x_i|x_j) using Markov transition rules forward and backward
%

% number of steps in chain
n=50;

% Set the transition rule!
p=0.05;
P=[(1-p) p ; 0 1];  

% initial state probabilities
pp0=0.99
pi1=[pp0; 1-pp0];
d=size(P,2);

% recursion for marginal probabilities
pmarg(1,1:d)=pi1';
for i=2:n,
    piSS=pmarg(i-1,:)*P;
    pmarg(i,1:d)=piSS;
end;    
figure(1); 
clf;
plot([1:n],pmarg(:,1),'k',[1:n],pmarg(:,2),'r--');
xlabel('Node index');
legend('p(x_i=1)','p(x_i=2)');
ylabel('Marginal probability');
axis([0 n+1 0 1]);

% Draw realizations 
real=10; % number of realizations
samp=zeros(n,real);
for j=1:real,
    vec=cumsum(pmarg(1,1:d)); % cumulative probability
    samp(1,j)=min(find(vec>rand)); % sample
    for s=2:n,
        vec=cumsum(P(samp(s-1,j),:)); % cumulative probability
        samp(s,j)=min(find(vec>rand)); % sample
    end;
end;    
sampV=samp;  

figure(2);
clf;
imagesc(sampV'); colormap(gray); 
%grid on; axis image;
xlabel('Node index');
ylabel('Realizations');


% VOI ANALYSIS 
Cost=10; % Cost
CostV=[0; 20]; % Cost if snow falls

% Post 
% Prior value
gg=max(-CostV'*pmarg',-Cost);
PV=sum(gg);
disp(sprintf('Prior value is %0.3g',PV));

% Evidence at single nodes , iobs=1....,n,
for iobs=1:n,
 
        % Low risk at iobs (last input is 1)
        Pcond=evidence_prop_Markov_chain(iobs,n,P,pmarg,1);
        if (iobs==20)
            figure(3);
            clf;
            subplot(2,1,1), title('High risk at 20');
            plot([1:n],Pcond(:,1),'k',[1:n],Pcond(:,2),'r--');
            xlabel('Node index');
            ylabel('Conditional probability');
        end;
        PoVobs1=max(-CostV'*Pcond',-Cost)*pmarg(iobs,1);
        % high risk at iobs (last input is 2)
        Pcond=evidence_prop_Markov_chain(iobs,n,P,pmarg,2);
        if (iobs==20)
           figure(3);
           subplot(2,1,2), title('High risk at 20');
           plot([1:n],Pcond(:,1),'k',[1:n],Pcond(:,2),'r--');
           xlabel('Node index');
           ylabel('Conditional probability');
        end
        PoVobs2=max(-CostV'*Pcond',-Cost)*pmarg(iobs,2);
        % VOI is expected value over data - prior value
        VOI(iobs)=sum(PoVobs1)+sum(PoVobs2)-PV;
end
    
figure(4);
clf;
plot(VOI,'k');
%st=title('VOI for one observation sites'); set(st,'FontSize',14);
ylabel('VOI');
xlabel('Test node j');
axis([1 n 0 1.1*max(VOI)]);


%% Subfile for evidence propagation

function Pcond=evidence_prop_Markov_chain(iobs,n,P,pmarg,xobs);


% Pcond is p(x_i=1|x_{iobs}=xobs) for i=1....n,
d=size(P,1);
Pcond=zeros(n,d);
if (xobs==1)
    % low risk at iobs
    Pcond(iobs,1)=1;
    pvecIOBS=[1;0];
else
    % high risk at iobs
    Pcond(iobs,2)=1;
    pvecIOBS=[0;1];
end;

% For higher indeces, forward recursion
pvec=pvecIOBS;
for i=iobs+1:n,
    pvec_n=pvec'*P;
    Pcond(i,1:d)=pvec_n;
    pvec=pvec_n';
end

% For smaller indeces, backward recursion
pvec=pvecIOBS;
for i=iobs-1:-1:1,
    % Backward transition matrix
    Pback=P'.*((1./pmarg(i+1,:)')*pmarg(i,:));
    pvec_n=pvec'*Pback;
    Pcond(i,1:d)=pvec_n;
    pvec=pvec_n';
end


