clear
echo on

%This script uses the attached function
%prob.m to calculate the probability of 
%exactly k occurrences out of n trials
%provided that the probability of one
%occurrence is p.

%Suppose we flip a coin 50 times.  What is
%the probability of getting heads 25 times?

x=prob(.5,25,50)

%which is actually pretty large.
pause
%We can also plot the probability distribution
%as a function of k:

figure(1)
k=[1:50];
plot(k,prob(.5,k,50))
xlabel('k')
ylabel('probability')
pause
%From the graph we see that the maximum
%probability occurs at k=25 as we might
%expect.
pause

%We can also look at the cumulative probability
%of getting k heads:

figure(2)
semilogy(k,cumsum(prob(.5,k,50)))
xlabel('k')
ylabel('probability of no more than k occurrences')
title('Coin Flip Probability Distribution')
grid on


%So the probability of getting no more than 10 heads
%or 10 tails out of 50 is about 10^-5, pretty
%small!

%Try playing around with this function to
%get a better feel for probability!
echo off


_________________________________________
cut here and store separately

function x=prob(p,k,n)

%This function calculates the probability of
%obtaining exactly k positive results out of
%N trials, given that the individualprobability
%of one positive outcome is given by p.  In this
%function p and n are scalars and k may be a
%vector.  The program is written so that overflow
%and underflow are avoided even for large values
%of n.  This is done by multiplying or dividing
%by some large number b to keep everything in 
%range.  The variable m keeps track of the number
%of times b has been used.
%
%call using x=prob(p,k,n)

b=10^10;
npts=length(k);

for i=1:npts
 xt=1;
 kt=k(i);
 m=0;

%This gets the first part of the binomial coef:
 for j=(n-kt+1):n
  xt=xt*j;
  if xt > b
   xt=xt/b;
   m=m+1;
  end
 end
 
%Now for the second part of the binomial coef:
 for j=1:kt
  xt=xt/j;
  if xt < 1/b
   xt=xt*b;
   m=m-1;
  end
 end
 
%And then multiplying by the probability of each
%individual occurrence:
 for j=1:kt
  xt=xt*p;
  if xt < 1/b
   xt=xt*b;
   m=m-1;
  end
 end
 
 for j=(kt+1):n
  xt=xt*(1-p);
  if xt < 1/b
   xt=xt*b;
   m=m-1;
  end
 end

%Finally we multiply by the quantity b^m to 
%correct for all of our adjustments:
 x(i)=xt*b^m;
end;