1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB Probability at Party function.

  1. Oct 8, 2013 #1
    1. The problem statement, all variables and given/known data
    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.


    2. Relevant equations



    3. The attempt at a solution
    Code (Text):

    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)
     
    Last edited by a moderator: Oct 9, 2013
  2. jcsd
  3. Oct 8, 2013 #2

    DrClaude

    User Avatar

    Staff: Mentor

    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##.
     
  4. Oct 8, 2013 #3
    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!
     
  5. Oct 9, 2013 #4

    DrClaude

    User Avatar

    Staff: Mentor

    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
    Code (Text):

    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.)
     
  6. Oct 9, 2013 #5
    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:
    Code (Text):
    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.
     
  7. Oct 9, 2013 #6

    DrClaude

    User Avatar

    Staff: Mentor

    I guess you mean q = 1.

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

    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.

    I imagine that this is a typo and that you have
    Code (Text):

      S = S * (1 - m/365) ;
      m = m+1;
     
     
  8. Oct 10, 2013 #7
    Got it! Apologies for the typos. Thanks for the help.
     
  9. Oct 10, 2013 #8
    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?
     
  10. Oct 10, 2013 #9

    DrClaude

    User Avatar

    Staff: Mentor

    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 =
    Code (Text):

      1.0000000
    + 0.0000001
    = 1.0000001
     
    whereas

    1.0000000 + 9 × 10-8 =
    Code (Text):

      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 × 105 + 10-7 =
    Code (Text):

      1.0000000 E 5
    + 0.000000000001 E 5
    = 1.0000000 E 5
     
    You have to consider instead (1 + ε) × 105, i.e., consider only the mantissa.
     
  11. Oct 10, 2013 #10
    Brilliant! Thanks for all the help!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: MATLAB Probability at Party function.
  1. MATLAB Function (Replies: 0)

Loading...