function treasure;

%
% Simple example of treasure 
%
%

% Probability success
pA=0.01;

% P positive test given distinction of interest is there
pBgivenA=0.95;
% P positive test given distinction of interest is NOT there
pBgivenAc=0.05;

% Bayes rule
%P(A|B)=P(A)P(B|A)/ ( P(A)P(B|A) + P(Ac) P(B|Ac) )
pAgivenB=(pA*pBgivenA)/(pA*pBgivenA+(1-pA)*pBgivenAc)


% Plot probability of rare event
pBgivenAcVec=[0.005:0.001:0.5];

for i=1:length(pBgivenAcVec),
    pBgivenAc=pBgivenAcVec(i);
    pAgivenBVec(i)=(pA*pBgivenA)/(pA*pBgivenA+(1-pA)*pBgivenAc);
end;

figure(10)
clf;
plot(pBgivenAcVec,pAgivenBVec,'k');
sy=ylabel('Conditional probability of success');
sx=xlabel('Probability of false positives');
%st=title('Rare event probabilities');
set(sx,'FontSize',16);
set(sy,'FontSize',16);
%set(st,'FontSize',16);
st=text(0.2,0.5,'p(x=1)=0.01, p(y=1|x=1)=0.95')
set(st,'FontSize',16);
stt=text(0.01,0.025,'p(x=1)=0.01');
set(stt,'FontSize',16);
st=text(0.1,0.4,'Plot shows p(x=1|y=1) for different p(y=1|x=0)')
set(st,'FontSize',16);
hold;
plot(pBgivenAcVec,0.01*ones(size(pBgivenAcVec)),'k--');
axis([0 max(pBgivenAcVec) 0 max(pAgivenBVec)])

p=pA;
Rev=100000;
Cost=[0:1000:100000];
gamma=[0.95 0.99 1];

%
% VOI computation
%
PV=zeros(size(Cost,2),1);
PostV=zeros(size(Cost,2),size(gamma,2));
VOI=zeros(size(Cost,2),size(gamma,2));
post0=zeros(size(Cost,2),size(gamma,2));
post1=zeros(size(Cost,2),size(gamma,2));
for i=1:size(Cost,2),
  PV(i,1)=max(0,Rev*p-Cost(i));
  for j=1:size(gamma,2),
    Pd0=(1-gamma(j))*p+gamma(j)*(1-p);
    Pd1=1-Pd0;
    post10=p*(1-gamma(j))/Pd0;
    post0(i,j)=max(Rev*post10-Cost(i),0);
    post11=p*gamma(j)/Pd1;
    post1(i,j)=max(Rev*post11-Cost(i),0);
    PostV(i,j)=Pd0*post0(i,j)+Pd1*post1(i,j);
    VOI(i,j)=max(PostV(i,j)-PV(i),0);
  end;
end;

figure(1)
clf;
plot(Cost,PV,'k'); 
      xlabel('Cost'); title('Prior value');
axis([0 2000 0 1000]);

figure(2)
clf;
plot(Cost,PV,'k'); 
      xlabel('Cost'); ylabel('value');
hold;
plot(Cost,PostV(:,3),'r');
%axis([0 2000 0 1000]);

figure(3)
clf;
plot(Cost,VOI(:,3),'k');
title('Value of information');
xlabel('Cost');

figure(4)
clf;
for j=1:size(gamma,2),
  subplot(1,size(gamma,2),j), plot(Cost,post1(:,j)','k'); hold; 
  plot(Cost,post0(:,j),'k-.'); plot(Cost,PV,'k--'); 
      xlabel('Cost'); title(sprintf('gamma is %0.5g',gamma(j)));
  ylabel('Values in different data cases and prior');  
end;

figure(4)
clf;
plot(Cost,VOI(:,1)','k'); hold; plot(Cost,VOI(:,2),'-.k');
      plot(Cost,VOI(:,3),'--k'); 
      sx=xlabel('Cost'); %st=title(sprintf('gamma is %0.5g',gamma(2)));
      %set(st,'FontSize',16);
  sy=ylabel('Value of information'); 
  set(sy,'FontSize',14);
  set(sx,'FontSize',14);
  sl=legend('VOI(y), \gamma=0.95','VOI(y), \gamma=0.99','VOI(x)');
  set(sl,'FontSize',16);
%axis([0 10000 0 120])


