% LMex030204_4th.m or LMex030308_3rd.m
% Exam 3.2.4, page 135-136 in 
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% and Larsen & Marx (2001) Third Edition Example 3.3.8 p. 139
% Written by E. Gallagher Eugene.Gallagher@umb.edu
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% 10/6/2010 revised 1/11/2010
% A major goal is to produce Figure 3.2.1, which will require a large
% number of points
% initialize arrays to store results
% Doomsday airlines has 2 aircraft --- a dilapidated 2-engine prob plane
% and an equally dilapidated under-maintained 4-enghine prop ploe. Each
% plane will land safely if at least half its engines are working properly.
% Given that you want to remain alive, under what conditions would you opt
% to fly in the two engine plane. Assume that the engines on each plane
% have the same probability of failing and that engine failures are
% independent events.
penginefails=zeros(999,1);
p2enginesafe=zeros(999,1);
p4enginesafe=zeros(999,1);
for i=1:9999
    penginefails(i)=i/10000;
    p2enginesafe(i)=sum(binopdf(1:2,2,1-penginefails(i)));
    p4enginesafe(i)=sum(binopdf(2:4,4,1-penginefails(i)));
end
plot(penginefails,p2enginesafe,'-r',...
    penginefails,p4enginesafe,'--g','LineWidth',3);figure(gcf)
xlabel('p=P(engine fails)','FontSize',16)
ylabel('P(Safe flight)','FontSize',16)
i=find(p2enginesafe>=p4enginesafe);
mi=min(i);
fprintf('Pick a 2 engine plane if P(enginefails) >= %8.6f \n',...
    penginefails(mi));
% A more precise estimate can be had using equation 3.2.3 in text, p 136
% and Matlab's built-in fzero function
Penginefails = fzero(@(p) (3*p-1)*(p-1), 0.001, 0.999);
fprintf('Now solve the problem more precisely using Matlab''s fzero.m\n')
fprintf('Pick a 2 engine plane if P(enginefails) >= %8.6f \n',...
    Penginefails);
P2enginesafe=sum(binopdf(1:2,2,1-Penginefails));
hold on
% This plot simply puts a vertical line in the graph
plot([Penginefails Penginefails],[0 P2enginesafe],'-.b','LineWidth',3);
legend('2-engine plane','4-engine plane','Location','NorthEast')
text(0.36,0.07,'1/3','FontSize',18)
title('Figure 3.2.1','FontSize',20)
figure(gcf)
hold off
