One thing not considered so far is that birthdays are not evenly distributed - see http://www.panix.com/~murphy/bday.html - and this affects the probability.
On Fri, Jan 27, 2012 at 9:50 AM, Linda Alvord <lindaalv...@verizon.net>wrote: > I'm glad that we have finally covered most of the conditions that could > occur. I'll enjoy pondering your work. I'm working on a "narrower" > challenge for next time that will not take us so far afield. > > Linda > > -----Original Message----- > From: programming-boun...@jsoftware.com > [mailto:programming-boun...@jsoftware.com] On Behalf Of Jose Mario > Quintana > Sent: Thursday, January 26, 2012 11:08 PM > To: Programming forum > Subject: Re: [Jprogramming] Challenge 4 Bountiful Birthdays > > Regarding the original task, one can proceed (knowing the actual densities) > via a shifted multinomial simulation, > > multinomial=. +/\ o [ <: o ((0: , [) I. ]) ? o $&0 o ] NB. dyadic verb > > 7 (samples=. ((densities o [) (2+multinomial) ]) ("0)) 10 10 10 10 10 NB. > day of the week samples > 2 3 4 6 5 6 2 6 2 6 > 3 4 2 8 2 5 4 5 3 5 > 3 2 2 6 2 3 4 5 2 3 > 4 5 5 6 2 6 5 7 2 3 > 5 2 3 5 5 5 4 5 4 7 > > mean=. +/ % # > > 365 (samplesmeans=. (mean"1 o samples)) 10000000 NB. day of the year 10 > million sample mean > 24.6141 > > 10 (] , mean) o (365 &samplesmeans) o # 500 NB. the original task > 24.264 23.782 24.334 24.516 25.016 24.704 25.05 24.514 23.93 25.25 24.536 > > (] , mean) o (365 &samplesmeans) o # f. NB. according to the "simple" > rules? > (] , +/ % #)@:(365&((+/ % #)"1@:((+/\^:_1@:(1 - */\@:(1 - ] %~ 1 + i.))@:[ > (2 + +/\@:[ <:@:((0: , [) I. ]) ?@:$&0@:]) ])"0)))@:# > > Regarding accuracy, among other things, it can be argued that the > distribution could even depend whether the experiment is conducted in the > northern or the southern hemisphere (see > http://www.panix.com/~murphy/bday.html and > http://answers.google.com/answers/threadview/id/280242.html). Models, > maps, > and other representations are ultimately doomed to be inaccurate; the > subject matter is not only too complex but also evolving; above all, of > course, my representation of "the world" that is my mind is affected as > well :) > > ________________________________________ > From: programming-boun...@jsoftware.com [programming-boun...@jsoftware.com > ] > On Behalf Of Jose Mario Quintana [josemarioquint...@2bestsystems.com] > Sent: Tuesday, January 24, 2012 2:23 PM > To: Programming forum > Subject: Re: [Jprogramming] Challenge 4 Bountiful Birthdays > > > +/(2+i.1000) * p * q NB. expected value > > 24.6166 > > I found the same solution in a slightly different way, > > ((2 + i.) +/ .* +/\^:_1@:((1 - */\)@:(1 - ] %~ 1 + i.))) 365 > 24.6166 > > The outline follows: > > It is easier to start dealing with the day of the week birthday process > first, > > (outcomes=. 2 + i.) 7 NB. all other outcomes have zero densities; thus, > they are irrelevant > 2 3 4 5 6 7 8 > > o=. @: > > (cp=. 1 - ] %~ 1 + i.) 7 NB. conditional probabilities the process will > not stop at each outcome given that it did not stop at its predecessor > 0.857143 0.714286 0.571429 0.428571 0.285714 0.142857 0 > > cdf=. 1 - */\ o cp NB. cumulative distribution function > > load'plot' > plot (0 0 , cdf) 7 NB. ploting the (smoothed) cdf > > densities=. +/\^:_1 o cdf NB. since cdf -: +/\ densities > > (mean=. outcomes +/ .* densities) 7 NB. formula for discrete densities > 4.01814 > > This generalizes to the day of the year birthday process, > > plot (0 0 , cdf) 365 > mean 365 > 24.6166 > > mean f. > (2 + i.) +/ .* +/\^:_1@:(1 - */\@:(1 - ] %~ 1 + i.)) > > > > ________________________________________ > From: programming-boun...@jsoftware.com [programming-boun...@jsoftware.com > ] > On Behalf Of Mike Day [mike_liz....@tiscali.co.uk] > Sent: Friday, January 20, 2012 7:54 PM > To: Programming forum > Subject: Re: [Jprogramming] Challenge 4 Bountiful Birthdays > > My "trial" function, listed earlier (and below) was > not quite correct, as it failed to count the > successful person. > > So it should be: > > trialb =: ([: # (] (,`]@.e.~ ([: ? 365"_)))^:_)"0 > > So we get, for example (but it's very slow! My variant > triala discussed with Linda is somewhat better): > > (mean, stdev) mean trialb 5000 100 $ _1 > 24.6133 0.180788 > > Linda thinks the mean should be somewhat lower, and > Brian thinks it's a lot lower. However, the standard > deviation suggests it's close to the true value. > > I think this is the way to find the true expected number > of people. We don't need Markov after all: > > Probability that (n-1) arrivals all have different > b/days: > > q =: Prod (1 - i%Y), 0<: i <: n-2, Y =~ 365 > > Probability that the nth arrival's b/day is one of > those present, ie is one of n-1 distinct bdays: > > p =: (n-1) % Y > > Expected value of number of arrivals for "success": > > Sum (2+i) pi * qi, 0 <: i <: n-2 > > In J: > 5{. q =: */\(1 - (365 %~(i.))) 1000 > 1 0.99726 0.991796 0.983644 0.972864 > > 5 {. p =: (365 %~>:@i. )1000 > 0.00273973 0.00547945 0.00821918 0.0109589 0.0136986 > > +/(2+i.1000) * p * q NB. expected value > 24.6166 > > This is not the same as the median, where the > probability q moves below 0.5, > > 21 22 { q > 0.524305 0.492703 > > As Roger observes, the index origin comes into play; > we should add one as the first person is 1, not zero (!) > and the median group size is therefore just below 23. > > This last is dealing with a slightly different problem: > what is the probability that a certain sized group of > people do (not) share a birthday? So we shouldn't be > surprised at the difference. > > Mike > > On 18/01/2012 3:17 PM, Mike Day wrote: > > People seem to be tackling two different problems. > > > > Variations on the Birthday Problem as I remember them: > > (a) what is the probability that two (or more) people > > share a birthday in a group of N people? > > (b) what should N be for the probability to be (say) 0.5 ? > > The somewhat counter-intuitive answers are dealt with in > > Roger's Wiki Essay, among many treatments, and also > > Pablo's message, below. The essential point is to > > consider the probability that there are no matches. > > > > However, Linda's single trial as stated is a random > > process with a stopping condition: > > take one person at a time until the new person shares a > > birthday with those already present. The result is the > > number of people including the new arrival. > > > > I expect you need a Markov Process approach to get the > > exact expected value for the stopping number. Not proved! > > > > Here's a stab at the required simulation, avoiding @ and @: > > though using [: > > > > NB. I use _1 as seed, so need to decrement the count > > > > trial =: (_1 + [: # (] (,`]@.e.~ ([: ? 365"_)))^:_)"0 > > > > trial 10#_1 NB. eg conduct 10 trials > > 27 19 29 2 24 42 30 9 34 33 > > > > mean =: +/%# NB. ok for vectors or columns of matrix > > > > ([:(;~mean) mean) TRIALS =: trial 500 10 $ _1 > > > > +-------+------------------------------------------------------------------- > + > > > > |23.5882|22.696 23.676 23.894 24.044 23.874 23.56 24.258 23.416 22.81 > > 23.654| > > > > +-------+------------------------------------------------------------------- > + > > > > > > These means are indeed close to N in problems > > (a) & (b) where the probability is ~0.50, namely > > 21 for 0.475695 and 22 for 0.507297, but not the > > same. > > > > I used 365 rather than Pablo's 365.25 . The simulation > > could be done for 365.25, using the integer 1461 (say). > > The stopping condition would be a bit more complicated. > > > > The deviation of trials is quite large: > > SS =: [: *: (-"1 mean) NB. squared deviations from mean > > stdev=: [: %: [: mean SS NB. Observed Standard deviation > > NB. not necessarily recommended for real, large sets of data > > > > (mean,:stdev) TRIALS > > 22.696 23.676 23.894 24.044 23.874 23.56 24.258 23.416 22.81 > > 23.654 > > 11.9378 12.6587 12.6917 12.5288 12.281 11.9741 12.1957 11.442 12.0969 > > 12.8718 > > > > NB. standard deviation of the means: > > > > (mean, stdev) mean TRIALS > > 23.5882 0.477041 > > > > Mike > > > > > > > > > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > -- Devon McCormick, CFA ^me^ at acm. org is my preferred e-mail ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
