%% The Bionomial Distribution

N=  100;        %Number of trials
p=  .02;        %Probability of success on each trial
q=(1-p);        %Probability of failure

%The proability mass function for k number of successes

P_k=inline('nchoosek(N,k)*p^k*(1-p)^(N-k)','k','p','N')

%Calculate P(k,p,N) for a range of k values
clear P
k_max=20;
k=(0:k_max);
for ii=1:length(k)
P(ii)=P_k(k(ii),p,N);
end

%Plot the probability distribution
figure
bar(k,P,'r','linewidth',2)
xlabel('Number of Successes (k)','fontsize',14)
ylabel('P(k,p,N)','fontsize',14)
s=sprintf('Binomial Distribution, p=%g, N=%d',p,N)
title(s)

%Check normalization
sum(P)

%Find the expectation of k (should be N*p)
k_exp=sum(k.*P)

%Find variance of k (should be N*p*(1-p)
k_var=sum( (k-N*p).^2.*P)

%% The Poisson Distribution

lambda=5;   %expected value

%The proability mass function for k number of successes
P_k=inline('lambda^k*exp(-lambda)/gamma(k+1)','k','lambda')

%Calculate P(k,p,N) for a range of k values
clear P
k_max=20;
k=(0:k_max);
for ii=1:length(k)
P(ii)=P_k(k(ii),lambda);
end

%Plot the probability distribution
figure
bar(k,P,'r','linewidth',2)
xlabel('Number of Successes (k)','fontsize',14)
ylabel('P(k,\lambda)','fontsize',14)
s=sprintf('Poisson Distribution, \\lambda=%g',lambda)
title(s)

%Check normalization
sum(P)

%Find the expectation of k (should be labmda)
k_exp=sum(k.*P)

%Find variance of k (should be labmda)
k_var=sum( (k-lambda).^2.*P)

%% The Gaussian (Normal) Distribution

sigma=  2;      %standard deviation
mu=     5;      %mean

%The proability density function 
P_x=inline('1/sqrt(2*pi*sigma^2)*exp( -(x-mu).^2/(2*sigma^2))','x','mu','sigma')

%Calculate P(x) for a range of x values
clear P
dx=.5;
x_max=20;
x=(-x_max:dx:x_max);
for ii=1:length(x)
P(ii)=P_x(x(ii),mu,sigma);
end

%Plot the probability density
figure
plot(x,P,'r','linewidth',2)
xlabel('x','fontsize',14)
ylabel('P(x)','fontsize',14)
s=sprintf('Normal Distribution, \\mu=%g, \\sigma=%g',mu,sigma)
title(s)

%Check normalization
sum(P*dx)

%Find the expectation of k (should be mu)
x_exp=sum(x.*P*dx)

%Find variance of k (should be sigma^2)
x_var=sum( (x-mu).^2.*P*dx)














