%% 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);