function [pom, pov]= BA_nornor2(data, sigma2, mu, tau2)
%-------------------------------------------------------------------
% [pom, pov]= BA_nornor(data, sigma2, mu, tau2)
% This function will plot the likelihood, prior, and the posterior
% in the Normal/Normal Case. Inputs are:
%    data - vector of observations
%    sigma2 - variance of data
%    mu, tau2 - prior mean/variance
% Output is posterior mean/variance  [pom, pov]
%    and the graphs of three densities.
% Example of use: 
%   > dat=10*randn(1,20);
%   > BA_nornor2(dat, 100, 20, 10)
%-----------------------------------------------------------------
n = length(data);
xbar = mean(data); 
set(0, 'DefaultAxesFontSize', 16);
fs = 16;
lw = 2.5;  
msize = 6;

figure(1)
xx = xbar-4*sqrt(sigma2):(sqrt(sigma2)/50):xbar+4*sqrt(sigma2);
yy = 1/(sqrt(2 * pi * sigma2/n)) .* exp(-1./(2 *sigma2/n).* (xx - xbar).^2 );
plot(xx,yy,'k-','linewidth',lw)
hold on
xxxx= mu -4*sqrt(tau2):(sqrt(tau2)/50):mu+4*sqrt(tau2);
yyyy = 1/(sqrt(2 * pi * tau2)) .* exp(-1./(2 *tau2).* (xxxx - mu).^2 );
plot(xxxx,yyyy,'k--','linewidth',lw)
hold on
xxx=mu - -4*sqrt(tau2):(sqrt(tau2)/50):mu+4*sqrt(tau2);
lambda = tau2/(tau2 + sigma2/n);
pom = lambda * xbar + (1-lambda)*mu;
pov = sigma2/n * tau2/(tau2 + sigma2/n);
xxx= pom -4*sqrt(pov):(sqrt(pov)/50):pom+4*sqrt(pov);
yyy = 1/(sqrt(2 * pi * pov)) .* exp(-1./(2 *pov).* (xxx - pom).^2 );
plot(xxx,yyy,'k:','linewidth',lw)
legend('likelihood','prior','posterior',1)
xlabel('\theta')
axis tight
hold off
print -depsc 'C:\NPBook\Bayes\bayesnorm.eps'