% from binomial to normal distribution
flag = 1;
switch flag
case 1
  n = 100; % number of events
  N = 10^6; % number of drawings
  p = 0.5;

  data = binornd(n,p,N,1);

  [heights,centers] = hist(data,30:1:70);
  heights = heights/sum(heights);
  
  mu = mean(data)
  sigma=sqrt(var(data))
  x = linspace(30,70);
  pdf = normpdf(x,mu,sigma);


  figure(1)
  plot(x,pdf,'r')
  hold on
  bar(centers, heights)
  hold off
  xlim([30,70])
  
  x = 0:n;
  pdf_Binomial = binopdf(x,n,p);
  pdf_Normal = normpdf(x,mu,sigma);
  % using exact values
  pdf_Normal = normpdf(x,n*p,sqrt(n*p*(1-p)));
  
  figure(2)
  plot(x,pdf_Binomial,x,pdf_Normal)
  legend('binomial' ,'normal')
  title('Binomial and normal distribution')
  figure(3)
  plot(x,pdf_Binomial-pdf_Normal)
  title('Difference of binomial and normal distribution')
  
case 2
  n = 10; % number of events
  N = 10^6; % number of drawings
  p = 0.2;
  data = binornd(n,p,N,1);

  [heights,centers] = hist(data,0:10);
  heights = heights/sum(heights); 
  mu = mean(data)
  sigma=sqrt(var(data))
  x = linspace(0,10);
  pdf = normpdf(x,mu,sigma);


  figure(1)
  plot(x,pdf,'r')
  hold on
  bar(centers, heights)
  hold off
  xlim([0,10])
  title('Histogram and normal distribution')
  
  x = 0:n;
  pdf_Binomial = binopdf(x,n,p);
  pdf_Normal = normpdf(x,mu,sigma);
  % using exact values
  pdf_Normal = normpdf(x,n*p,sqrt(n*p*(1-p)));
  
  figure(2)
  plot(x,pdf_Binomial,x,pdf_Normal)
  legend('binomial' ,'normal')
  title('Binomial and normal distribution')
  figure(3)
  plot(x,pdf_Binomial-pdf_Normal)
  title('Difference of binomial and normal distribution')
  end%case