%%%clt

clear all; clf reset;

np=4;  %values of the r.v. are 1,..., np

for (n=1:60)

  pause(0.2);clf;
  p=zeros(1,np);

  %uniform
  for (i=1:np)
    p(i)=1/np;
  end
  %custom
  p(1)=0.5; p(2)=0.1; p(3)=0.2; p(4)=0.2;

  %p(1)=0.5; p(2)=0; p(3)=0.1; p(4)=0.4;

  %p(1)=0.8; p(2)=0; p(3)=0.05; p(4)=0.15;

  ex=[1:np]*p';
  ex2=([1:np].^2)*p';
  sig=sqrt(ex2-ex^2);

  %convolution
  p1=p;
  np1=np;

  for (i=2:n)
    np2=np1+np;
    p2=zeros(1,np2);
    for (j0=1:np)
      for (j1=1:np1)
        p2(j0+j1)=p2(j0+j1)+p(j0)*p1(j1);
      end
    end
    p1=p2; np1=np2;
  end

  hold on;
  bar([1:np1], p1);

  %plot normal density with same exp. and var.
  mp=n*ex;
  sp=sig*sqrt(n);
  lb=floor(n*ex-5*sig*sqrt(n));
  ub=floor(n*ex+5*sig*sqrt(n));

  inp=[lb:0.01:ub];

  xlim([lb, ub]);
  yub=1.1*max(p1);
  ylim([0,yub]);
  %xticks([lb:2:ub]);
  xlabel('$k$','FontSize',15, 'interpreter', 'latex');
  ylabel('$P(S_n=k$)','FontSize',15, 'interpreter', 'latex');
  str =sprintf('n=%d',n);
  text(lb+1,yub,str,'FontSize',20);

  plot(inp,  exp(-(inp-mp).^2/(2*sp^2))/(sqrt(2*pi)*sp), 'LineWidth', 2, 'Color', [1,0,0]);
  if (n<=4) waitforbuttonpress; end
  hold off;
end
