MATLAB Probability at Party function.

[SherlockOhms]

Homework Statement

So, you're asked to write a function that computes the probability of two people at a party having the same birthday. The probability function itself is

P(n) = 0, n = 1

P(n) = 1 - (capital pi from k = 1 to k = n - 1)(1 - k/365), 2 <=n <= 365

P(n) = 1, n >= 366

The function is used to calculate the minimum value of m (your chosen n) for which P(m) >= q. Where q (0, 1], i.e. an arbitrary probability entered by the user. The function should accept q as an input and return m as an output.

Homework Equations

The Attempt at a Solution

function [m] = Prob(q)
m = 0;                                  
P = 0;         
                              
while P < q                             
    m = m + 1;                       
    if m == 1                          
        P = 0;

    elseif m >= 2 && m <= 365         
        S = 1;                               
        for k = 1: m - 1 
            S = S * (1 - k/365);
        end 
        P = 1 - S;   
                                            
    else                        
        P = 1;  
    end 
end

Essentially, I want the program to accept a value of q and the to run until the value of P is greater than or equal to q. Once the the value of P is returned, m will have it's minimum value for which P >= q due to m being incremented during each iteration of the while loop. However, for some reason the value of "m"s being returned is too low and I can't seem to see why that is. For example, when I enter a value of 1 into the function, m is returned as 153 an not 366. Can anyone spot why this may be? I've looked through the code a few times now and can't seem to spot why. Is the incrementation being skipped or is there a problem with the nested for loop which calculates the value of P that involves the product of 1 - k/365 from k = 1, to k = m - 1? Thanks for any possible comments in advance. (Sorry, had to use underscores for indenting)

 

[DrClaude]

Plot the value of S(m), i.e.,
$$
S(m) = \prod_{k=1}^{m-1} \left( 1 - \frac{k}{365} \right)
$$
as a function of m on a logarithmic plot, and then look up machine epsilon.

Basically, $P = 1 - S(m) \approx 1$ for $m \ge 153$.

 

[SherlockOhms]

Right. So, it's basically just rounding P to 1 at m = 153, so the while loop's condition is met. Is there anyway in which I can alter the condition in the while loop to avoid this rounding problem? Thanks!

 

[DrClaude]

MATLAB won't work on my computer just now, so I can't test this, but my guess is to transform the test from
$$
P(m) = 1- S(m) < q
$$
to
$$
S(m) > 1 - q
$$
which would look something like

function [m] = Prob(q)
  m = 1; 
  
  if q > 0
    S = 1

    while S > 1-q && m <= 365  
      m = m + 1; 
      
      S = S * (1 - (m - 1)/365);
    end
  end

end

(Note: please use the CODE tags when posting code. It preserves the spacing and makes it easier to copy and paste.)

 

[SherlockOhms]

Thanks for that. So, I guess I could test for q being equal to 0 that I could then say m = 366.
Also, I have an alternate program written out that actually works correctly it goes like this:

function m = Prob (q)

if (q <=0) || (q >1)
error ('Input q out of bounds : q must lie in (0 ,1]. ')
end

if q==1
 m =366;

else
 m = 1;
 S = 1;
 while 1- S < q
  S = S * (1-S /365) ;
  m = m+1;
 end
end

Why is it that 1 - S is not rounded to 1, causing the condition to be met early, in this case as S → eps? Thanks.

 

[DrClaude]

Thanks for that. So, I guess I could test for q being equal to 0 that I could then say m = 366.

I guess you mean q = 1.

function m = Prob (q)

if (q <=0) || (q >1)
error ('Input q out of bounds : q must lie in (0 ,1]. ')
end

I didn't mention it previously, but checking parameters is a very good (you could say necessary) practice.

 while 1- S < q

Why is it that 1 - S is not rounded to 1, causing the condition to be met early, in this case as S → eps?

This will give you the same problems of rounding due to the machine epsilon. If S gets too small, then 1-S = 0 and the condition is met. But you have put in an extra condition, testing for q = 1. So the biggest value of q≠1 you can input is $q = 1 - \epsilon$, such that you are comparing $1-S < q$ up to $S=\epsilon$, hence it works.

  S = S * (1-S /365) ;
  m = m+1;

I imagine that this is a typo and that you have

  S = S * (1 - m/365) ;
  m = m+1;

 

[SherlockOhms]

Got it! Apologies for the typos. Thanks for the help.

 

[SherlockOhms]

Could you explain the role of eps? Is it that it's the smallest gap possible between two floating point numbers therefore it's impossible to perform arithmetic with any number less than eps? Is this in anyway correct or have I just convinced myself of something totally wrong?

 

[DrClaude]

Could you explain the role of eps? Is it that it's the smallest gap possible between two floating point numbers therefore it's impossible to perform arithmetic with any number less than eps? Is this in anyway correct or have I just convinced myself of something totally wrong?

That sounds correct. To understand exactly how it works in a computer, we would have to look are floating points numbers are handled, but you would get the same with a calculator. If your calculator can display 8 digits (and no rounding is ever performed), then ε = 10^-7

1.0000000 + 10^-7 =

  1.0000000
+ 0.0000001
= 1.0000001

whereas

1.0000000 + 9 × 10^-8 =

  1.0000000
+ 0.00000009
= 1.0000000

When considering calculations in a program, the value of ε is of relative, in the sense that it will vary depending on the exponent,
1 × 10⁵ + 10^-7 =

  1.0000000 E 5
+ 0.000000000001 E 5
= 1.0000000 E 5

You have to consider instead (1 + ε) × 10⁵, i.e., consider only the mantissa.

 

[SherlockOhms]

Brilliant! Thanks for all the help!