Monte Carlo simulation for probability of a multi part machine working

% each A part works 98% of the time. Each B works 91% of the time.
% The machine has 12A and 7B parts. To function, at least 5B and 12A 
% parts must function. What is the probability to function. 
% The simulation does not work, it just gives me errors or outputs 0 no
% matter how I phrase the syntax. What is wrong with the simulation?
N = 1000000 ;
Parts= zeros(N,1);
RelativefreqWork=zeros(N,1);
countA=0;
countB=0;
Aworks=0;
Bworks=0;
m=0;
for i=1:1000000
A=randi(12*100,1); %makes 12 random numbers between 1-100
B=randi(7*100,1); %makes 7 random numbers between 1-100
if A(i)<=98; %probability to see if random number A is <98
    countA=countA+1; %counts every time A is <98
    if countA>=12 %if all 12 A parts work
        Aworks=Aworks+1 %count +1 into machine section A working
    end
if B(i)<=91; %probability to see if random number B is <91
    countB=countB+1; %counts every time number B is <91
    if countB>=5 %if at least 5 B parts function
        Bworks=Bworks+1 %counts +1 into machine section B working
    end
        if Bworks>=5 && Aworks>=12 %if both conditions are met
           m=m+1;
        end
end
RelativefreqWork(i)=m/i
end
end

%% Plot Histogram
figure
sturgesn=ceil(log2(length(m))+1);
[x,c]=hist(m,sturges); %prints the hisogram
h=gca;
set(h,'XTick',c);
x/sum(x);
bar(c,x/sum(x));
xlabel('Probability that machine works'), ylabel('Relative Frequency');

%% convergence plot
figure;
semilogx(1:N,(1-(m)));
grid on
axis tight;
xlabel('Number of trials'),ylabel('Prob that it is mean');

Answer 1

[Alan Stevens]

Here's one possibility:

% each A part works 98% of the time. Each B works 91% of the time.
% The machine has 12A and 7B parts. To function, at least 5B and 12A 
% parts must function. What is the probability to function. 
% The simulation does not work, it just gives me errors or outputs 0 no
% matter how I phrase the syntax. What is wrong with the simulation?
N = 10^5 ;

RelativefreqWork=zeros(N,1);
Working = zeros(N,1);

for i=1:N
    A = rand(12,1);           % 12 numbers between 0 and 1
    B = rand(7,1);            %  7 numbers between 0 and 1
    iA = find(A>0.02);        % list of working A's
    iB = find(B>0.09);        % list of working B's         
    Aworks = numel(iA);       % number of working A's
    Bworks = numel(iB);       % number of working B's
    if Aworks==12 && Bworks>=5
       Working(i) = 1;        
    end
    RelativefreqWork(i) = sum(Working)/i;
end

disp(RelativefreqWork(N))

    0.7689

%% convergence plot
figure;
plot(1:N,RelativefreqWork);
grid on
axis([0 N 0.5 1]);
xlabel('Number of trials'),ylabel('Prob machine is working');

Answer 2

[Yongjian Feng]

Hello Petch,

First of all, check the help of randi. You can type help randi from the command line window of matlab. You didn't use it correctly.

Secondly, you need to loop through your 12 simulation of A and 7 simulation of B. You had only one loop here for 1000000.

Thanks,

Yongjian

Answer 3

[Image Analyst]

This is the way I'd approach it. I adapted your code:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;

% Each A part works 98% of the time. Each B works 91% of the time.
% The machine has 12 A and 7 B parts. 
% User said "To function, at least 5 B and 12 A parts must function."
% In other words: To function, at least 12 A parts (all of them) and 5 B parts must function. 
% What is the probability to function?
numMachines = 1000000; % A million
countA = zeros(numMachines, 1);
countB = zeros(numMachines, 1);
countBoth = zeros(numMachines, 1);
for k = 1 : numMachines % For evary machine...
	% Construct parts A and B for this particular machine.
	A = rand(12, 1); %makes 12 random numbers between 1-100
	B = rand(7, 1); %makes 7 random numbers between 1-100
	
	% See which of the 12 A parts work.
	workingPartsA = A <= 0.98; % Probability to see if random number A is < 98%
	% See which of the 7 B parts work.
	workingPartsB = B <= 0.91; % Probability to see if random number B is < 91%
	
	% See if the number of A working parts is 12 because ALL must work.
	if sum(workingPartsA) == 12 
		countA(k) = 1;
	end
	% See if the number of B working parts is at least 5 out of the 7.
	if sum(workingPartsB) >= 5
		countB(k) = 1;
	end
	% See if both work
	if countA(k) && countB(k)
		countBoth(k) = 1;
	end
end
% Print out how many worked.
fprintf('Out of %d machines, %7d (%7.4f%%) worked because both part A and part B worked.\n',...
	numMachines, sum(countBoth), 100 * sum(countBoth)/numMachines);
% Print out how many failed due to A alone failing.
numAAlone = sum(~countA & countB); % A didn't work but B did work.
fprintf('Out of %d machines, %7d (%7.4f%%) worked because part A alone failed.\n',...
	numMachines, numAAlone, 100 * numAAlone/numMachines);
% Print out how many failed due to B alone failing.
numBAlone = sum(countA & ~countB); % A worked but B didn't.
fprintf('Out of %d machines, %7d (%7.4f%%) worked because part B alone failed.\n',...
	numMachines, numBAlone, 100 * numBAlone/numMachines);
% Print out how many failed due to both A and B failing.
numBothFailed = sum(~countA & ~countB); % A and B both failed.
fprintf('Out of %d machines, %7d (%7.4f%%) worked because both parts A and B failed.\n',...
	numMachines, numBothFailed, 100 * numBothFailed/numMachines);

You get:

Out of 1000000 machines, 769438 (76.9438%) worked because both part A and part B worked.

Out of 1000000 machines, 211277 (21.1277%) worked because part A alone failed.

Out of 1000000 machines, 15158 ( 1.5158%) worked because part B alone failed.

Out of 1000000 machines, 4127 ( 0.4127%) worked because both parts A and B failed.

You can plot bar charts to show the relative amounts if you want.