%% Bayes estimate, see book PCV Figur 4.1
close all
clearvars
% 1D Gaussian density
g=@(x,m,s) exp(-1/2*(x-m).^2/s^2)/sqrt(2*pi*s^2);
% uniform density
u = @(x,a,b) 1/(b-a)*(x>a & x < b);
% chosen discretization of parameter theta
dx = 0.01;
% range of investigation
range = -3:dx:6;


%% set variables of demo
% prior, smaller mean 
mu1 = 1.2
si1 = 0.1
% larger mean
mu2 = 1.6
si2 = 0.1
% observational noise
sil = 0.15
% asymmetry of prior
P1 = 0.55;
% range of uniform density of outliers
r_min = 0;
r_max = 6;

for k=1:2
    if k==1 % inlier and non-robust likelihood
        P_out = 0;
        l=1.5;
        % likelihood = Gaussian
        plx=g(range,l,sil);
        r_ind = 350:600;
        % mixture prior
        px=P1*g(range,mu1,si1)+(1-P1)*g(range,mu2,si2);
        % posterior, scaled
        pxl = px.*plx;
        pxl = pxl/sum(pxl)/dx;
    else
        l=2.2; % outlier and robust likelihood
        P_out = 0.2;
        % likelihood = mixture of Gaussian and uniform
        plx=(1-P_out)*g(range,l,sil)+P_out*u(range,r_min,r_max);
        % mixture prior
        px=P1*g(range,mu1,si1)+(1-P1)*g(range,mu2,si2);
        % posterior, scaled
        pxl = px.*plx;
        pxl = pxl/sum(pxl)/dx;
    end
    % boundary for plot
    xl_max = max([mu1+4*si1,mu2+4*si2,l+4*sil]);
    xl_min = max([mu1-4*si1,mu2-4*si2,l-4*sil]);
    
    % find maximum of posterior
    [~,indx]=max(pxl);
    xmax=range(indx)
    
    
    if k==1
        figure
        plot(range(r_ind),px(r_ind),'-k','LineWidth',2)
        xticks(0.6:0.2:2.6)
        title(['p(\theta)'])
        xlim([0.5,xl_max])
        ylim([0,5])
        
        figure
        plot(range(r_ind),plx(r_ind),'--k','LineWidth',2)
        xticks(0.6:0.2:2.6)
        title('p(y|\theta)')
        xlim([0.5,xl_max])
        ylim([0,5])
        
        figure
        plot(range(r_ind),pxl(r_ind),'-k','LineWidth',1)
        xticks(0.6:0.2:2.6)
        title(['p(\theta|y)'])
        xlim([0.5,xl_max])
        ylim([0,5])
        
        xlim([0.6,xl_max])
        ylim([0,5])
        
        figure
        hold on
        plot(range(r_ind),px(r_ind),'-k','LineWidth',2)
        plot(range(r_ind),plx(r_ind),'--k','LineWidth',2)
        plot(range(r_ind),pxl(r_ind),'-k','LineWidth',1)
        xticks(0.6:0.2:2.6)
        %title(['Thick: p(\theta), dashed: p(y|\theta), thin: p(\theta|y). y=',num2str(l),',\theta=',num2str(xmax)])
        
        title(['Thick: p(\theta), dashed: p(y|\theta), thin: p(\theta|y).  y=',num2str(l),',\theta=',num2str(xmax)])
        xlim([0.6,xl_max])
        ylim([0,5])
        
    else
        figure
        plot(range,px,'-k','LineWidth',2)
        title(['p(\theta)'])
        xlim([0.5,3.5])
        ylim([0,5])
        
        figure
        hold on
        plot(range,plx,'--k','LineWidth',2)
        plot([-1,0],[0,0],'--k','LineWidth',2)
        plot([6,7],[0,0],'--k','LineWidth',2)
        title(['P(out)=',num2str(P_out),'%,p(y|\theta,possible-outlier)'])
        xlim([-1,7])
        ylim([0,5])
        
        figure
        plot(range,pxl,'-k','LineWidth',1)
        title(['P(out)=',num2str(P_out),'%,p(\theta|y)'])
        xlim([0.5,3.])
        ylim([0,5])
        
        xlim([0.5,3.5])
        ylim([0,5])
        
        figure
        hold on
        plot(range,px,'-k','LineWidth',2)
        plot(range,plx,'--k','LineWidth',2)
        plot(range,pxl,'-k','LineWidth',1)
        title(['P(out)=',num2str(P_out),...
            '%. Thick: p(\theta), dashed: p(y|\theta), thin: p(\theta|y). y=',num2str(l),',\theta=',num2str(xmax)])
        xlim([0.5,3.5])
        ylim([0,5])
    end
    
end