% Lecture 3A example: Study the distribution of a sum
%
% Parameters:
%    n       Number of terms in the sum
%    minx    Min value of each term (default=1)
%    maxx    Max value of each term (default=6)
%    px      Probabilities of the term values (default = uniform)
%    vis     Boolean; if true, a bar plot is shown
%    r       Width of the plot (r standard deviations left and right),
%            default is to show full scale.
%
% Try minx=0, maxx=1 for indicator variables.
%
% Returns:
%    ps      The probabilities of sum values
%    s       The possible sum values
%    Es      Expexted value of the sum
%    Ds      Standard deviation of the sum
% 
function [ps,s,Es,Ds]=L3Asumdist(n, minx, maxx, px, vis, r)
if nargin<1, n=2; end
if nargin<2, minx=1; end
if nargin<3, maxx=6; end
if nargin<4, k=maxx-minx+1; px=repmat(1/k, 1, k); end
if nargin<5, vis=false; end
if nargin<6, r=nan; end

% Start from one term
ps = px;

% Apply convolution formula n-1 times
for k=2:n
    ps = conv(ps, px);
end

s = (n*minx):(n*maxx);
Es = sum(s.*ps);
Ds = sqrt(sum((s-Es).^2 .* ps));

if vis
    bar(s,ps, 'facecolor', 'b')
    xlabel('s')
    ylabel('P(S=s)')
    title(sprintf('%d terms',n))
    grid on

    if not(isnan(r))
        % Scale picture to show only 4 standard deviations
        % left and right from the mean.
        set(gca,'xlim',[Es-r*Ds Es+r*Ds]);
    end
end