Probability Function of people at a party.

I've recently asked a question about this particular section/piece of code. I've made a bit of preogress but have ran into some trouble. I'll restate the question and code in question. 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. Here's my code:

function [m] = Prob(q)
m = 0; %assign m initial value of 0
P = 0; %assign P a value of 0 for now so that while loop will start
while P < q %runs until P >= q
      m = m + 1; %increments m as long as condition above is true
      if m == 1 %checks first condition
          P = 0;
      elseif m >= 2 && m <= 365 %checks second condition
          S = 1; %variable that stores the product of (1 - k/365) from 1 -> m-1
          for k = 1: m - 1 
              S = S * (1 - k/365);
          end
          P = 1 - S; %calculates P by subtracting S(the total product) from 1                                      
      else %checks third condition
        P = 1; %assigned 1 if m >= 365
      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.

Answer 1

[Walter Roberson]

Suppose S < eps(1), then P = 1 - S would be P = 1 due to roundoff.

Try testing 1-q against S instead of 1-S against q

[David]

Gotcha. That would make sense seeing as at about the k = 153 mark, I'd assume 1 - k/365 will have gotten quite small. As in, test 1 - q < S for the condition in the while loop? Thanks. Also, would S not be rounded to 0 in the above condition if S < eps?