% L5Bdirichlet   Demo of Dirichlet distribution (for three alternatives)
%
% By default we have observed counts [5 3 2].
% By default we use 100 * 100 grid.
%
% Hint: Try bigger values, e.g. [50 30 20] for a 100-person sample
% or [500 300 200] for a 1000-person sample, typical for opinion polls.
%
% With large samples the credible region will be much smaller.
% With high counts you may need to increase the grid resolution (N).
%
function L5Bdirichlet(counts, N)

if nargin<1, counts = [5 3 2]; end
if nargin<2, N      = 100; end

clf

grid = linspace(0,1,N);
[p,q] = meshgrid(grid,grid);
r = 1-p-q;
r(r<0) = 0;
P = [p(:) q(:) r(:)];

% Prior
subplot(1,3,1)
A = [1 1 1];    % parameters for uniform Dirichlet distribution
fprior = dirpdf(P, A);
fprior = reshape(fprior, N, N);
%surf(p, q, fprior);
surf(p,q,fprior, 'linewidth', 0.2)
colorbar
title('Prior density')
xlabel('p')
ylabel('q')

% Posterior
subplot(1,3,2)
B = 1+counts;
fpost = dirpdf(P, B);
fpost = reshape(fpost, N, N);
%contour3(p, q, fpost);
contour3(p,q,fpost,15,'k','linewidth',2);
hold on;
surf(p,q,fpost, 'Edgecolor', 'none')
colorbar
title('Posterior density');
xlabel('p')
ylabel('q')


% Posterior 50%, 95% and 99% credible regions
subplot(1,3,3)
I = credibleregion(fpost, 0.50);
J = credibleregion(fpost, 0.95);
K = credibleregion(fpost, 0.99);
imagesc(grid, grid, I+J+K + (fprior>0))
set(gca,'ydir','normal')
title('50%/95%/99% credible regions');
xlabel('p')
ylabel('q')