function [x,y] = poissonsim(mean, tries) % % POISSON distribution generator % function [x,y] = gen(mean, tries) % % mean = the mean of the POISSON distribution % tries = number of times the data was taken % x = the vector of simulated Poisson distribution -- which will be % different everytime this generator is run % y = the binned simulated Poisson distribution % % This function plots errorbars only if mean <=10. Otherwise, plot % gets to crowded. If you want to edit it, copy it and poisson.m to % your home directory and change it. % first need to find the probability integral... gran = 1; % this is the granularity of the distribution % realistically set to 1 m=0; j = 0; pset = 0:m; pint = 0:m; x = 0:j; randpick = 0:m; if (mean > 4) && (mean < 25) maxb = 4*ceil(mean); else if mean >= 25 maxb = 2.5*ceil(mean); else maxb = 16; end end for i = 0:gran:maxb m = m + 1; pset(m) = poisson(i,mean); pint(m) = sum(pset(1:m))*gran; end i = 0:gran:maxb; bin = ones(maxb/gran+2,1)-1; % generate set of totally random numbers for j = 1:tries randpick(j) = rand; % now need to find which bins those random numbers correspond to tf = 0; k = 1; while tf == 0 if randpick(j) < pint(k) bin(k) = bin(k) + 1; tf = 1; x(j) = k-1; else k = k + 1; % this bin is overflow if k > maxb/gran+1 bin(k) = bin(k) + 1; tf = 1; x(j) = k-1; end end end end y = bin; % now plot di = 0:gran/10:maxb; %set(axes,'fontsize',14); plot(di,poisson(di,mean)*tries*gran,'b-', 'linewidth',2) hold on plot(i,bin(1:length(bin)-1),'r.','markersize',20) err = sqrt(bin(1:length(bin)-1)); if mean <=10 errorbar(i,bin(1:length(bin)-1),err,'r.') end xlabel('Frequency bins','fontsize',16) ylabel('Number of events','fontsize',16) title('Poisson simulator','fontsize',18) axis([0 maxb 0 1.5*tries*poisson(mean,mean)]) hold off