coins.m

% collected data (1 = definitely fair, 7 = definitely unfair)
street = [2; 3; 5; 2; 4; 5; 2; 6; 7];
streetMum = [2; 3; 6; 2; 1; 7; 2; 2; 7];
bank = [1;1;4;1;3;5;1;3;6];
bankMum = [2;2;6;2;3;6;2;2;6];
labels = {'Coin 1' 'Coin 2' 'Coin 3' 'Coin 4' 'Coin 5' 'Coin 6' 'Coin 7' 'Coin 8' 'Coin 9'};

% number of heads and tails for each coin flip sequence
coins = num2cell([3 2; 1 4; 5 0; 4 6; 8 2; 0 10; 10 16; 21 5; 26 0], 2);

% manually trained logistic parameters
a = 0.8;
b = -2.5;

% run prediction on each flip sequence
predictions = cellfun(@(x) predict(a, b, x(1), x(2)), coins);
% scale to the 1-7 scale used for collected data
% (1 = fair, 7 = unfair)
scaled = arrayfun(@(x) 1+6*(1-x), predictions);

% plot
axis = 1:9;
plot(axis, scaled, axis, street, axis, streetMum, axis, bank, axis, bankMum)
set(gca, 'XTick', axis, 'XTickLabel', labels)
legend('Model', 'Street (me)', 'Street (friend)', 'Bank (me)', 'Bank (friend)')

logposterior.m

% computes the log posterior odds ratio
% log P(H_1|D)/P(H_2|D) =
%   log P(D|H_1)/P(D|H_2) + log P(H_1)/P(H_2)
% we assume priors are equal, so the last term on the
% RHS simplifies to 0
function l = logposterior(H, T)
    l = log(PDH1(H, T)/PDH2(H, T));
end

PDH1.m

% computes P(D|H_1)
% H_1: the fair coin hypothesis
% heads and tails are equally likely, so P = (1/2)^n where n is coin flip
% sequence length (= H + T, where H is the number of heads and T the number
% of tails)
function P = PDH1(H, T)
    P = 1/(2^(H+T));
end

PDH2.m

% computes the likelihood P(D|H_2)
% H_2: the weighted coin hypothesis
function P = PDH2(H, T)
    syms n
    % sum P(D|theta_n) * P(theta_n|H2)
    % (thetan ranges uniformly from 0 to 1)
    % this approximates \int_0^1 P(D|theta) P(theta|H_2) dtheta
    Pf = symsum(PDTheta(thetan(n), H, T)*PThetaH2(thetan(n)), n, 0, 100);
    P = double(Pf);
end

function theta = thetan(n)
    theta = 1/100*n;
end

% computes P(Theta | H_2)
% H_2: the weighted coin hypothesis (for all possible theta)
% (theta is the weighting parameter for the coin)
function P = PThetaH2(theta)
    P = 1/100;
end

% computes P(D | Theta) for a specific weighted coin
% (theta is the weighting parameter for the coin)
% (H is the number of heads, T is the number of tails) 
function P = PDTheta(theta, H, T)
    P = PHTheta(theta)^H * PTTheta(theta)^T;
end

% computes P(Heads | Theta) for a specific weighted coin
% (theta is the weighting parameter)
function P = PHTheta(theta)
    P = theta;
end

% computes P(Tails | Theta) for a specific weighted coin
% (theta is the weighting parameter)
function P = PTTheta(theta)
    P = 1 - theta;
end

predict.m

% pass the log posterior odds ratio through a logistic
% function with parameters a, b
function f = predict(a, b, H, T)
    f = 1/(1 + exp(-a*logposterior(H, T)+b));
end