%
% Solution to exercise 1, homework 1, part 2
% CAAM 415, 10/16/08
%

%number of release sites
n1 = 200;
x = 0:1:n1; %used to compute pdfs and plot them

p1 = [0.01 0.05 0.1];

figure('Position',[200 200 1000 500]);
for i = 1:3
  %simply go through all the cases one by one and compute
  l1 = p1(i)*n1; %quantal size
  by1 = binopdf(x,n1,p1(i)); %compute binomial 
  py1 = poisspdf(x,l1); %compute poisson
  bpy = [by1;py1]';
  
  subplot(1,3,i);
  h = bar(x,bpy); %plot everything
  set(gca,'XLim', [-1 40]);
  xlabel('number of quanta');
  ylabel('probability of release');
  legend(h,'binomial', 'Poisson');
  
  dbpy = diff(bpy,1,2);
  %subselect the binomial values larger than .05
  binds = find(by1 >= 0.05); 
  %maximal relative error
  max_err_rel1(i) = max(abs(dbpy(binds))./by1(binds)'); 
end;


info_str = sprintf('N = 200, relative errors: %.3f %.3f %.3f',max_err_rel1);
disp(info_str);

n2 = 4;
x = 0:1:n2;


p2 = [0.1 0.5 0.75];

figure('Position',[200 400 1000 500]);
for i = 1:3
  %simply go through all the cases one by one and compute
  l2 = p2(i)*n2; %quantal size
  by2 = binopdf(x,n2,p2(i)); %compute binomial 
  py2 = poisspdf(x,l2); %compute poisson
  bpy = [by2;py2]';
  
  subplot(1,3,i);
  h = bar(x,bpy); %plot everything
  set(gca,'XLim', [-1 5]);
  xlabel('number of quanta');
  ylabel('probability of release');
  legend(h,'binomial', 'Poisson');
  
  dbpy = diff(bpy,1,2);
  %subselect the binomial values larger than .05
  binds = find(by2 >= 0.05); 
  %maximal relative error
  max_err_rel2(i) = max(abs(dbpy(binds))./by2(binds)'); 
end;


info_str = sprintf('N = 4, relative errors: %.3f %.3f %.3f',max_err_rel2);
disp(info_str);