%% A modern Introduction to Probability and Statistics
% Chapter 2 - Outcomes, events and probability.
%
% Special thanks to David Dobor (http://david-dobor.github.io/2033/) great 
% examples used for this lecture.

%% Event: Toss a fair coin.
% Probability of tossing a fair coin is 0.5 for both head and tail.
% P(H) = P(T) = 1/2.
% 1 - Heads, 0 - Tails

x = rand <= 0.5 % command > has an 'if-like' structure, it return logical
% 0 or 1 when condition is satisfied.

%% Event: Roll a die
% Probability of tossing a fair die is 1/6 for each side of the die.
% P(1) = P(2) = P(3) = P(4) = P(5) = P(6) = 1/6.

x = ceil(6*rand)

%% Event: Roll a pair of dice
x = ceil(6*rand) + ceil(6*rand)

%% Toss a coin a bunch of times
n = 10
x = rand(1, n) <= .5

%% Toss a coin a bunch of times; count the number of heads
% 10 tosses of a fair coin. 1 - Heads, 0 - Tails

n = 10
t = rand(1, n) <= .5
x = sum(t)

%% Toss a coin a bunch of times in many trials
clear
p = .7
n = 100
trials = 10000
x = zeros(1,trials);

for i = 1:trials
    t = rand(1,n) <= p;
    x(i) = sum(t);
end


hist(x, 0:n)
title('Distribution of dice rolls');
xlabel('Counts of Heads');
ylabel('Trials');

%% Event: The birthday experiment
% Let's ask 30 people randomly what is the month they were born in?
%
% First we should model this event by distributing a probability for each
% month by how many days it has (assuming everyone were not born on a leap
% year).

daysInYear = 365;
months = {'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', ...
    'Aug', 'Sep', 'Oct', 'Nov', 'Dec'}
% Another (neat) way of generating vector of month names:
monthsIndex = 1:1:24
months = datestr(datenum(2000,monthsIndex,1),'mmm')
daysInMonth = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];

% Calculate probability of a person being born in each month!
birthMonthProbabilities = daysInMonth/daysInYear;

% Lets plot the probability
figure
bar(birthMonthProbabilities)

axis([0 13 0 inf])
title('Probability to be born on a month');
xlabel('Months');
ylabel('Probabilities');
set(gca, 'XTick', 1:12)
set(gca, 'XTickLabel', months);

% And now lets see what happened when we asked a few people on which month
% are they born
numberOfPeople = 30;
peopleAnswered = [];

for person = 1:numberOfPeople
    aux = ceil(rand * 365);
    if aux >= 1 && aux < 31
        bornIn = 1;
    else if aux >= 31 && aux < 59
            bornIn = 2;
        else if aux >= 60 && aux < 90
                bornIn = 3;
            else if aux >= 91 && aux < 120
                    bornIn = 4;
                else if aux >= 121 && aux < 151
                        bornIn = 5;
                    else if aux >= 152 && aux < 181
                            bornIn = 6;
                        else if aux >= 182 && aux < 212
                                bornIn = 7;
                            else if aux >= 213 && aux < 242
                                    bornIn = 8;
                                else if aux >= 243 && aux < 273
                                        bornIn = 9;
                                    else if aux >= 274 && aux < 304
                                            bornIn = 10;
                                        else if aux >= 305 && aux < 334
                                                bornIn = 11;
                                            else if aux >= 335 && aux < 365
                                                    bornIn = 12;
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end



%    bornIn =  round(unifrnd(1,12,1,1))
peopleAnswered = [peopleAnswered, bornIn];
end

peopleAnswered

% Now let's see when people are usually born
figure
hist(peopleAnswered)

axis([0 13 0 inf])
title('Distribution of months people are born in');
xlabel('Months');
ylabel('Probabilities');
set(gca, 'XTick', 1:12)
set(gca, 'XTickLabel', months);

% And lastly let's compare the two
figure
subplot(1,2,1) % sublot allows us to plot multiple plots on same figure.
bar(birthMonthProbabilities)
axis([0 13 0 inf])
title('Probability to be born on a month');
xlabel('Months');
ylabel('Probabilities');
set(gca, 'XTick', 1:12)
set(gca, 'XTickLabel', months);
subplot(1,2,2)
hist(peopleAnswered)
axis([0 13 0 inf])
title('Distribution of months people are born in');
xlabel('Months');
ylabel('Probabilities');
set(gca, 'XTick', 1:12)
set(gca, 'XTickLabel', months);

%% Event: Deal a Poker Hand
deck = randperm(52);
hand = deck(1:5)

%% Event: Probablities of Poker Hands

% model counts of possible outcomes for each hand in Poker
% Hint: http://en.wikipedia.org/wiki/Poker_probability
straightFlush=40;
fourOfaKind=13*48;
fullHouse=13*12*4*nchoosek(4,2);
flush=4*nchoosek(13,5)-40;
straight=10*4^5-40;
threeOfaKind=13*4*48*44/2;
twoPair=nchoosek(13,2)*nchoosek(4,2)*nchoosek(4,2)*44;
pair=13*nchoosek(4,2)*48*44*40/factorial(3);
squat=nchoosek(13,5)*4^5-straight-flush-straightflush;

hands=[straightFlush, fourOfaKind, fullHouse, flush, straight, ...
    threeOfaKind,twoPair,pair,squat];
% sum of all Poker hands outcomes used to get probability for each hand
total=sum(hands)
% check if the total = totalshouldbe
totalShouldBe = nchoosek(52,5)
format long
probabilities=hands/total
format short

%Lets plot the probabilities
handsName = {'Straight flush','Four of a kind','Full house','Flush', ...
    'Straight','Three of a kind','Two pair','Pair','Squat'};
figure
bar(probabilities)
axis([0 length(handsName)+1 0 inf])
title('Probability to get a poker hand');
xlabel('Poker Hands');
ylabel('Probabilities');
set(gca, 'XTick', 1:length(handsName)+1)
set(gca, 'XTickLabel', handsName);