calculate probability in MATLAB

THIS IS A MULTINOMIAL PROBABILITY DISTRIBUTION PROBLEM!

say we are given 4 outcome below along with the probabilities

Calculate the probability such that

my attempted code is

% Define the probabilities for each outcome
HH = 0.013;
HL = 0.267;
LH = 0.031;
LL = 0.689;

% Initialize the total probability
totalProbability = 0; 

% Loop through possible counts for HL and LH
for nHL = 0:20
    for nLH = 0:2
        nHH = 0; 
        nLL = 100 - nHH - nHL - nLH

       probability = mnpdf([nHH  nHL  nLH nLL], [HH   HL  LH  LL])
       totalProbability = totalProbability + probability;
        end
end

% Display the result
disp(['The total probability is: ', num2str(totalProbability)])

However the code gives the value 0.0056902 but the real answer is 0.00132705 so clearly something is wrong with my code.

[Torsten]

Since

% Define the probabilities for each outcome
HH = 0.013;
HL = 0.267;
LH = 0.031;
LL = 0.689;
% Initialize the total probability
totalProbability = 0; 
% Loop through possible counts for HL and LH
for nHH = 0:100
  for nHL = 0:100
    for nLH = 0:100 
       nLL = 100 - nHH - nHL - nLH;
       if nLL >= 0
         probability = mnpdf([nHH  nHL  nLH nLL], [HH   HL  LH  LL]);
         totalProbability = totalProbability + probability;
       end
    end
  end
end
format long
totalProbability

totalProbability = 
   0.999999999999955

% Display the result
disp(['The total probability is: ', num2str(totalProbability)])

The total probability is: 1

I'd say that your code is correct. Where did you get 0.00132705 as result from ?

[N/A]

https://eng.libretexts.org/Bookshelves/Industrial_and_Systems_Engineering/Chemical_Process_Dynamics_and_Controls_(Woolf)/13%3A_Statistics_and_Probability_Background/13.10%3A_Multinomial_Distributions

[Torsten]

The argument in MATHEMATICA's CDF - function on the web page for the example given is wrong.

What is calculated on the web page is P(n1=0,n2<=20,n3<=2,n4<=78) which equals P(n1=0,n2=20,n3=2,n4=78) (because 0 + 20 + 2 + 78 is the only combination such that the sum equals 100):

format long
HH = 0.013;
HL = 0.267;
LH = 0.031;
LL = 0.689;
p = [HH HL LH LL];
x = [0 20 2 78];
mnpdf(x,p)

ans = 
   0.001327047052481

But what is needed is P(n1=0,n2<=20,n3<=2,n4<=100).

I have no access to MATHEMATICA, but my guess is that the correct command should be

<< Needs["MultivariateStatistics`"]
<< multinomial = MultinomialDistribution[100, {0.013, 0.267, 0.031, 0.689}]
<< CDF[multinomial, {0, 20, 2, 100}]

Answer 1

[Nipun]

Hi Nafisa,

I understand that you are working on a multinomial probability distribution problem and need help with your MATLAB code. Based on the shared information, I recommend using the following looping and assignment of frequencies.

% Define the probabilities for each outcome
HH = 0.013;
HL = 0.267;
LH = 0.031;
LL = 0.689;

% Initialize the total probability
totalProbability = 0;

% Loop through possible counts for HL and LH
for nHL = 0:20
    for nLH = 0:2
        nHH = 0;
        nLL = 100 - nHH - nHL - nLH;
        
        % Check if nLL is within valid bounds
        if nLL >= 0
            probability = mnpdf([nHH nHL nLH nLL], [HH HL LH LL]);
            totalProbability = totalProbability + probability;
        end
    end
end

% Display the result
disp(['The total probability is: ', num2str(totalProbability)])

This code ensures that the frequency nLL is non-negative and correctly sums up the probabilities for all valid combinations of nHL and nLH.

Hope this helps.

Regards,

Nipun