% ECE 248, UCSB, Homework 1 P1.7 a) and b) SOLUTION (due Oct. 22, 2010)
% Katie Byl
fprintf('ECE 248, UCSB, Homework 1 P1.7 a) and b) SOLUTION. Fall 2010.\n')
clear all; format compact % clear everything... print compactly

pwin = .55;   % probability of winning a single bet
BRstart = 2;  % Start with a bankroll ("BR") of 2 betting units
Bet = 1;      % Dollar value of a single bet
nbets = 1000;  % number of bets to simulation
Prob_BR = zeros(nbets+1,nbets+BRstart+1);  % bins for probability mass distribution
Prob_BR(1,BRstart + 1) = 1;  % 100% probability we have this money.

for n=1:nbets
    for bi = 1:(BRstart+n)  % bi: "bet index"
        bi_lose = bi-1;  % losing reduces index by one
        bi_win = bi+1;   % winning increases index by one
        if bi == 1  % however, bi=1 means we are already bankrupt!
            bi_win = 1;  % goto the min bin: bi=1 (bankrupt)
            bi_lose = 1; % goto the min bin: bi=1 (bankrupt)
        end
        Prob_BR(n+1,bi_win) = Prob_BR(n+1,bi_win) + ...
            pwin*Prob_BR(n,bi); % new bin if you WIN
        Prob_BR(n+1,bi_lose) = Prob_BR(n+1,bi_lose) + ...
            (1-pwin)*Prob_BR(n,bi); % new bin if you LOSE
    end
end
Money_each_bin = [0:(nbets+BRstart)]';
cum_EV = Bet*Prob_BR*Money_each_bin;
Xlist = [0:length(Prob_BR(:,1))-1];

% **------------* PLOT PROBABILITY OF RUIN vs N BETS *----------------**
figure(1); clf; subplot(211);  
plot(Xlist,Prob_BR(:,1),'.-'); grid on; axis([0 nbets 0 1])
title('HW1p7-a-i : Probability of RUIN')
xlabel('After the nth bet');  ylabel('Probability you are now BANKRUPT')
p_ruin_str = num2str(Prob_BR(end,1),'%.12f');
fprintf('\nPart A)\nAfter %d trials, the estimated probability of ruin is: %s\n',...
    nbets, p_ruin_str);
text(nbets,Prob_BR(end,1)+.1,['$P_{ruin} \approx $', p_ruin_str, '~~~~'],...
    'interpreter','latex','HorizontalAlignment','right','FontSize',18)
% **----------* PLOT EXPECTED VALUE OF BANKROLL vs N BETS *-----------**
subplot(212); 
plot([0 nbets],BRstart*Bet+(pwin-(1-pwin))*Bet*[0 nbets],'b--'); hold on
plot(Xlist',cum_EV,'r.-'); grid on
title('HW1p7-a-ii : Cumulative Winnings: expected value')
xlabel('After the nth bet'); ylabel('Expected winnings (EV), in betting units')
legend('If bankrupcy cannot occur','Actual bankroll-limited case',...
    'Location','NorthWest')
% ********* PART B: 
% Pbr(2) = Pbr(1)*Pbr(1) (Probability of going bankrupt from 1 vs 2
% bets as your initial bankroll.)  So we can express Pbr(1) as:
%     Pbr(1) = Plose*1 + Pwin*(Pbr(1)*Pbr(1))
%     0 = pwin*P^2 - P + (1-pwin)  [rewriting the algebra above]
R = roots([pwin, -1, (1-pwin)]);   % Solve the algebra above...
P_ruin_partB = min(R);             % Take the value not equal to 1
fprintf(['\nPart B)\nP(1) = %.11f is the probability of ruin, given you\n', ...
    'start with ONE betting unit.\n'],P_ruin_partB)
fprintf(['P(2) = P(1)^2 = %.11f is the probability of ruin given you start\n',...
    'with 2 units. This is (as expected) close to the value in part A.\n\n'],...
    P_ruin_partB^2)
    
% ===== END OF FILE =====