function Problem2;

% 
% Code for parts of Problem 2:
%
% Univariate Gaussian example. 
% pdf
% p(x)=N(mu,sigma^2)=(1/sqrt(2*pi*sigma^2))*exp(-(x-mu)^2/(2*sigma^2))
%
% x is the profit of project. 
% 
% File computes the VOI of perfect information. 
% Comparison of analytical and  Monte Carlo solutions.
%

% PROPERTIES:
% (x-mu) is a straight line which is negative x<mu and positive for x>mu
% The Gaussian pdf is symmetric around mu. This means
%\int (x-mu) p(x) dx = 
% \int_{-infty}^{mu} (x-mu) p(x) dx +\int_{mu}^{infty} (x-mu) p(x) dx =0
%
% Standard normal has mu=0, sigma=1:
% phi(z)=(1/sqrt(2*pi))*exp(-0.5*z^2). 
% dphi /dz = - z * phi(z)
% 
% phi(z)=phi(-z).  (normpdf(z) )
% Cumulative function Phi(z)=1-Phi(-z). (normcdf(z) )
%


% Parameter settings
% Mean
mu_ref=0;
mu=[-2:0.1:2];
% Standard deviation
sigma_ref=1;
sigma=[0.1:0.1:2];

% Number of Monte Carlo samples
B=1000;

%
% sigma_ref, varying mean
% ------------------------

% Analytical, Value of perfect information
VOI_anal=mu.*normcdf(mu./sigma_ref)+sigma_ref*normpdf(mu./sigma_ref)-max(mu,0);

% Monte Carlo samples
xb=ones(B,1)*mu+sigma_ref*randn(B,1)*ones(1,length(mu)); % xb=normrnd(mu,sigma,B,1);
% Monte Carlo approximation of VOI(x)
VOI_MC=mean(max(xb,0),1)-max(mean(xb,1),0);

figure(1);
clf;
plot(mu,VOI_anal),
hold;
plot(mu,VOI_MC,'k');
legend('Analytical','Monte Carlo');
xlabel('prior mean project'); ylabel(' VOI ');
title(sprintf('VOI versus prior mean. Reference sigma is %0.3g',sigma_ref));

%
% mu_ref, varying standard deviation
% ------------------------

% Analytical, Value of perfect information
VOI_anal=mu_ref*normcdf(mu_ref./sigma)+sigma.*normpdf(mu_ref./sigma)-max(mu_ref,0);

% Monte Carlo samples
xb=ones(B,1)*mu_ref*ones(1,length(sigma))+randn(B,1)*sigma; % xb=normrnd(mu,sigma,B,1);
VOI_MC=mean(max(xb,0),1)-max(mean(xb,1),0);

figure(2);
clf;
plot(sigma,VOI_anal),
hold;
plot(sigma,VOI_MC,'k');
legend('Analytical','Monte Carlo');
xlabel('standard deviation project'); ylabel(' VOI ');
title(sprintf('VOI versus prior uncertainty. Reference mean is %0.3g',mu_ref));


