function Monday_project;

%
% Two projects are Gaussian distributed a priori. 
% Computing VOI when
% - two IMPERFECT test (both projects, Gaussian likelihood)
% - one IMPERFECT test (one of the projects, Gaussian likelihood)
% 
% CALLS SUBFILE AT END:
% voi_calc_GaussProj(mu,rho,tau,F);

% p(x)=N(0,Sigma), Sigma=[1 rho;rho 1]
% p(y|x)=N(x,T), T=tau^2*I_m, m is 1 (partial) or 2 (total). 
%
% p(x|x) \propto p(x)p(y|x)
% \propto exp( -0.5*x'*(Sigma\x)-0.5*x'*(T\x)+x'*(T\y) )
% = exp (-0.5*x'* ( Sigma^{-1}+T^{-1}) * x +x'*(T\y))   )
%
% Conditional recognized by matching with
% -0.5 * x'* Sigma^{-1}_{x|y} x + x'*Sigma^{-1}_{x|y} mu_{x|y}
%
% p(y)=N(E(E(y|x)), Var(E(y|x))+E(Var(y|x)))=N(0,Sigma+T)
%

% vector of correlations
rhoV=[0.01:0.01:0.99];

tau=1; % measurement error / high accuracy

% Prior mean is set to 0.
mu=zeros(2,1);

% Partial test
F=[1 0]; % Picks first project.


% Initialize VOI for total and partial test
VOI_1=zeros(size(rhoV));
VOI_p1=zeros(size(rhoV));
% loop over correlations
for i=1:length(rhoV),
    
    % Correlation
    rho=rhoV(i);

    % TOTAL TEST 
    VOI_1(i)=voi_calc_GaussProj(mu,rho,tau);
    
    % PARTIAL TEST
    VOI_p1(i)=voi_calc_GaussProj(mu,rho,tau,F);

end
  
% Plotting VOI versus correlation
figure(1)
clf;
plot(rhoV,VOI_1,'k.',rhoV,VOI_p1,'b--');
hold;
legend('Total test','Partial test','Location','SouthEast');
ylabel('VOI');
xlabel('Correlation parameter');


% Decision regions
% argmax ( VOI2-(P1+P2), VOI-P1, None)
% Boundaries: P2=VOI2-VOI1, P2=VOI2-P1, P1=VOI1
%

% set correlation 
rho=0.7; % project correlation 

% TOTAL TEST
VOI2=voi_calc_GaussProj(mu,rho,tau);

% PARTIAL TEST
VOI1=voi_calc_GaussProj(mu,rho,tau,F);

VOI1

VOI2
   
% price range
figure(2);
clf;
plot([VOI1 VOI1],[VOI2-VOI1 2*(VOI2-VOI1)],'k');
hold;
plot([0 VOI1],[VOI2-VOI1 VOI2-VOI1],'k');
plot([VOI1 VOI2],[VOI2-VOI1 0],'k');
text(VOI1/3,(VOI2-VOI1)/3,'Total');
text(VOI1/2,1.1*(VOI2-VOI1),'Partial');
text(VOI1*1.1,(VOI2-VOI1),'No testing');
xlabel('Price of first test');
ylabel('Price of second test');

figure(3);
clf;
plot([VOI1 VOI1],[VOI2-VOI1 2*(VOI2-VOI1)],'k');
hold;
plot([0 VOI1],[VOI2-VOI1 VOI2-VOI1],'k');
plot([VOI1 VOI2],[VOI2-VOI1 0],'k');
text(VOI1/3,(VOI2-VOI1)/3,'Total');
text(VOI1/2,1.1*(VOI2-VOI1),'Partial');
text(VOI1*1.1,(VOI2-VOI1),'No testing');
xlabel('Price of first test');
ylabel('Price of second test');
% Color plotting to check lines
P1=[VOI1/20:VOI1/10:2*VOI1];
P2=[VOI1/20:VOI1/10:10*(VOI2-VOI1)];
for i=1:length(P1),
    for j=1:length(P2),
        if (VOI2<(P1(i)+P2(j)))
            if (VOI1<P1(i))
                plot(P1(i),P2(j),'ko');
            else
                plot(P1(i),P2(j),'rx');
            end
        else
            if (VOI1-P1(i)>VOI2-P1(i)-P2(j))
                plot(P1(i),P2(j),'rx');
            else
                plot(P1(i),P2(j),'b*');
            end
        end
    end
end


% ----Subfile: analytical VOI. Gaussian model. Free selection  --------

function VOI=voi_calc_GaussProj(mu,rho,tau,F);


if (nargin==3)
    F=eye(2);
end;

Sigma=[1 rho;rho 1];
T=tau^2*eye(size(F,1));
S=Sigma*F'*((F*Sigma*F'+T)\(F*Sigma));  
Sv=sqrt(diag(S));
if (Sv(2)>0)

    VOI=mu(1)*normcdf(mu(1)/Sv(1))+Sv(1)*normpdf(mu(1)/Sv(1))+...
        mu(2)*normcdf(mu(2)/Sv(2))+Sv(2)*normpdf(mu(2)/Sv(2))-sum(max(mu,0));

    % Note, when mu=zeros(2,1) entries cancel
else
    
    % No correlation, rho=0, no variance reduction for second project.
    % No value of information for second project
    % Avoid division on 0
    VOI=mu(1)*normcdf(mu(1)/Sv(1))+Sv(1)*normpdf(mu(1)/Sv(1))+sum(max(mu(1),0));
   %
end;


% Note : Matrix identity:
% 
% Scond = (Sigma^{-1}+T^{-1})^{-1}
%      = Sigma*(Sigma+T)^{-1}*T
%      = T*(Sigma+T)^{-1}*Sigma...
%      =  Sigma-Sigma*(Sigma+T)^{-1}*Sigma,
%
% This means that :
%
%    Scond=inv(inv(Sigma)+inv(T));
%    Salt=(Scond/T)*(Sigma+T)*(T\Scond)
%    S=Sigma*((Sigma+T)\Sigma);
%    S and Salt are identical
%
%  With partial information F=[1 0]:
%
%    Scond=inv(inv(Sigma)+F'*inv(T)*F);
%    Salt=((Scond*F')/T)*(F*Sigma*F'+T)*(T\(F*Scond))
%    S=Sigma*F'*((F*Sigma*F'+T)\(F*Sigma));
%    S and Salt are identical
%
% Proving conditional results is often easier with Scond. 
% VOI calculations are oftern easier with S rather than Salt.
%
