%  Using MATLAB to work with probability density functions  
%  for introducing pdf's in Brazil course
%    students don't have Statistics Toolbox
%    for  use with MATLAB 7 
%    with TOOLBOXES
%     SYMBOLIC MATH
%     STATISTICS
%   modified May 2008
%   by Carol Horvitz

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB lessons:
%  1) disttool
%  2) pdf 
%  3) various built-ins
%      a) normpdf
%      b) exppdf
%      c) poisspdf
%      d) binopdf
%      e) unidpdf
%  
%  4) writing your own non-uniform discrete probability function
%  7) if ... else ... end  loops

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all

help pdf
% we can choose one of interest and get the equation


% let's try the normal distribution with the example from Case, p.434
%which parameters do you need?
help normpdf

x=-4:.1:8; %set the range of outcomes for the evaluating the function
mu=2.5                   % set the mean
sigma=(2.25)^(0.5) % set the standard deviation

y1=pdf('norm',x,mu,sigma) % one way to get it
y2=normpdf(x,mu,sigma) %another way to get it

figure(1)
plot(x,y1); hold on;   %plot this line and keep it
plot(x,y2,'ro');            %add another line, indicate this line by red open circles
hold off
title(['Normal distribution, \mu =' num2str(mu) ', \sigma = ' num2str(sigma)])
xlabel('outcome')
ylabel('probability density')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% let's try the exponential distribution with the example from Case, p.434
%which parameters do you need to set?
clear all
close all

help exppdf

x=0:.1:3;                  %set the range of outcomes for the evaluating the function
mu=1/2.0                   % set the mean

y1=pdf('exp',x,mu) % one way to get it
y2=exppdf(x,mu) %another way to get it

figure(2)
plot(x,y1); hold on;   %plot this line and keep it
plot(x,y2,'ro');            %add another line, indicate this line by red open circles
hold off
title(['Exponential distribution, \mu =' num2str(mu) ', \lambda = ' num2str(1/mu)])
xlabel('outcome')
ylabel('probability density')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% let's try the lognormal distribution
%The normal and lognormal distributions are closely related. 
%If X is distributed lognormal with parameters µ and sigma, 
% then ln(X) is distributed normal with parameters µ and sigma. 

%The lognormal distribution is applicable 
%  when the quantity of interest must be positive,
%  since ln(X) exists only when the random variable X is positive. 
%Economists often model the distribution of income 
%   using a lognormal distribution.
%Cumulative population size is distributed this way.
%   (when something multiplicative causes the distribution...)

%which parameters do you need to set?
clear all
close all

help lognpdf %for more details


mu=40000;                     % set the mean
sigma=1.25;                    % set the standard deviation
x = (10:1000:125010)';      % set the range of x

figure(20)
y = lognpdf(x,log(mu),sigma);
plot(x,y);
text(mu,.1E-5,['\downarrowmean = ' num2str(mu)],'FontSize',16)
set(gca,'xtick',[0 30000 60000 90000 120000])
set(gca,'xticklabel',str2mat('0','30,000','60,000',...
                             '90,000','120,000'))
title(['Lognormal Distribution, \mu =' num2str(mu) ', \sigma = ' num2str(sigma)])
xlabel('outcome')
ylabel('probability density')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% let's try the Poisson distribution with the example from Case, p.433
%which parameters do you need to set?
clear all
close all

help poisspdf

x=0:10;                    %set the range of outcomes for the evaluating the function
lambda=4                  % set the mean (which is equal to the variance)

y1=pdf('poiss',x,lambda) % one way to get it
y2=poisspdf(x,lambda) %another way to get it

figure(3)
bar(x,y1); hold on;    %plot this as bars, since it is discrete and keep the bars
plot(x,y2,'ro');            %add these dots, indicate by red open circles
hold off
title(['Poisson distribution, \lambda = ' num2str(lambda)])
xlabel('outcome = number of successes')
ylabel('probability')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% let's try the Binomial distribution with the example from Case, p.433
%which parameters do you need to set?
clear all
close all

help binopdf

x=0:10;                    %set the range of outcomes for the evaluating the function
n=10                        % set the number of trials
p= 0.1                      % set the probability of "success" for each trial

y1=pdf('bino',x,n,p) % one way to get it
y2=binopdf(x,n,p) %another way to get it

figure(4)
bar(x,y1); hold on;    %plot this as bars, since it is discrete and keep the bars
plot(x,y2,'ro');            %add these dots, indicate by red open circles
hold off
title(['Binomial distribution, number of trials = ' num2str(n) ' , p =' num2str(p)])
xlabel('outcome = number of successes')
ylabel('probability')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% let's try a Uniform discrete distribution with an example, six-sided unweighte die
% which parameters do you need to set?
clear all
close all

help unidpdf

x=1:6;                       %set the range of outcomes for the evaluating the function
n=size(x,2);                %how many there are determines the probability of each
y1=pdf('unid',x,n); % one way to get it
y2=unidpdf(x,n); %another way to get it
for i=1:n             %another way to get it
y3(i)=1/n;
end

figure(5)
bar(x,y1); hold on;    %plot this as bars, since it is discrete and keep the bars
plot(x,y2,'ro');            %add these dots, indicate by red open circles
plot(x,y3,'g+');           % add these plus signs, indicate by green circles
hold off
title(['Uniform discrete distribution, no. of possible outcomes = ' num2str(n) ])
xlabel('outcome')
ylabel('probability')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%what about discrete, non-uniform probability a weighted six sided die?
clear all
close all

x=1:6;                       %set the range of outcomes for the evaluating the function
n=size(x,2);                %how many there are determines the probability of each

w= 2                                  %weighting
p_weighted=w/n                 %probability of weighted side
p_regular=[1-(w/n)]/[n-1]   %probability of the other sides
weighted_choice = 2            %let  2 be the weigted side

% find the probability for each outcome 
for i=1:n
    if i==weighted_choice
y1(i)=p_weighted;   %weight the die
    else
        y1(i)=p_regular;
    end
end

y1_cumprob=cumsum(y1)

figure(5)
subplot(2,1,1)
bar(x,y1);    %plot this as bars, since it is discrete and keep the bars
title(['Non-uniform discrete distribution, weighted side is ' num2str(weighted_choice) ])
xlabel('outcome')
ylabel('probability')
subplot(2,1,2)
bar(x,y1_cumprob)
xlabel('given outcome + all outcomes to the left')
ylabel('cumulative probability')

%we will use this type of construct in the "choosestudent" script,
% but first:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Consider the uniform continuous distribution
clear all
close all

a=2    %the lowest possible outcome
b=6    %the highest possible outcome

%what is the probability density function for x

% if x < a, the pdf=0
% if x > b, the pdf=0

% if x is between a and b, inclusive...what is its value?

% Draw a picture to help you answer this question...
%    the shape is a rectangle

% Its area = 1, its width =b-a,, its area= width*height
%   its height must be

fvector=0:0.01:10                                      %set some x values to look at
inrange = ((fvector >= a) & (fvector <= b)); %use logical operators
outofrange = ((fvector < a) | (fvector > b));
y1(inrange) = (1/(b-a)); %assign probabilities accordingly
y1(outofrange) = 0;


figure(6)
plot(fvector,y1,'r:','LineWidth', 2');    %plot this as bars, since it is discrete and keep the bars
title(['Uniform continuous distribution, lowest possible outcome= ' num2str(a) ', highest possible outcome =' num2str(b) ])
xlabel('outcome')
ylabel('probability density')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% another thing to do is
%call up the interactive pdf tool in the Statistics Tool Box
disttool
%use it to investigate a few distributions
% for each one, notice how many parameters need to be set
% we will answer some questions in class... take notes on how to use this
% interactive tool

% now for "chooseastudent"....