▼
Posts: 33
Threads: 12
Joined: Jan 1970
I'm trying to establish a geometric progression formula for the following series but I'm not sure it's a geometric progression. I need to be able to solve the nth term of the following series:
(364/365) * (363/365) * (362/365)...
I want to use 'a' for the changing numerator and 'b' for unchanging denominator and r for the ratio. I know the formula for a geometric progression but it's not that simple. What is the answer here? I need a formula that will fit all similar conditions that I can program on an HP 41. Right now I'm getting the answer by looping but I think there's a better way. Thanks.
▼
Posts: 536
Threads: 56
Joined: Jul 2005
how about (aPr)/b^r for r terms. where aPr means permutations of a taking r. ?
Edited: 15 June 2004, 6:51 p.m.
Posts: 97
Threads: 25
Joined: Jan 1970
Hello,
try this:
FACTORIAL(n-1)/n^(n-1) = GAMMA(n)*n^(1-n)
in your case n=365
Csaba
▼
Posts: 97
Threads: 25
Joined: Jan 1970
FACTORIAL(n-1)/n^(n-1) = GAMMA(n)*n^(1-n)
I'm so sorry, I think my mind is not workin today... And I can't to delete my message... Sorry again!
Cs.
▼
Posts: 33
Threads: 12
Joined: Jan 1970
The problem with factorials is that calculators and computers can't handle super large factoials. They overflow. Besides I'm looking for a partial or a truncated factorial. Logs of factorials may be able to handle that. With your proposed solution, when a = 364, b = 365 and n = 23, do you get .493? If not the real problem isn't solved. Try the problem manually, doesn't take that long, and verify the .493 answer and then formulate it for any value of a & b & n, and then you will feel comfortable with working on a solution.
Posts: 2,448
Threads: 90
Joined: Jul 2005
Quote:
(364/365) * (363/365) * (362/365)...
Won't this series go on for 364 cycles, decreasing all the time, and at an increasingly rapid rate, until you get to (0/365) at whcih point it goes disjointly to zero?
You could very easily write a loop, using i as a counter variable, and loop 364 times.....what am I missing here?
Regards,Bill
▼
Posts: 33
Threads: 12
Joined: Jan 1970
"You could very easily write a loop, using i as a counter variable, and loop 364 times.....what am I missing here?"
You're under the impression that I haven't solved my problem. Oh, I've solved it. I have a loop routine that works very well on a 41 but I can't loop with a 19BII and so I want to use some kind of progression formula that will find partial factorials. At n = 23 the ratio becomes .493. I'm beginning to think there is no formula. I haven't heard back from Dr. Math yet.
▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Jim:
Jim posted:
"I can't loop with a 19BII and so I want to use some kind of progression formula that will find partial factorials. At n = 23 the ratio becomes
.493. I'm beginning to think there is no formula."
Yes, there is. A "partial" factorial is technically called a Pochhammer Symbol, which itself is nothing more than a Gamma function divided over another gamma function, like this:
(x)n = x*(x+1)*...*(x+n-1) = Gamma(x+n)/Gamma(x)
This is the clue for solving your problem non-iteratively.
In HP-71B BASIC (i.e.: nearly plain-vanilla English), your problem would be solved like this:
10 B=365 @ N=23
20 DISP FACT(B)/FACT(B-N)/B^N
where FACT is the factorial function. However, upon running this results in an overflow, as FACT(365) exceeds the range for real numbers in the HP-71B and in most any other calculator or system.
What can we do, then ? Easy. Just use the asymptotic series for the factorial function, aka Stirling's series, to compute the natural logarithm of x! instead of x!. This computation won't overflow for your arguments, and once we have the logarithm, we just use the exponential function EXP (aka e^x) to deliver the final result. Stirling series is:
x! = (x/e)^x*Sqrt(2*Pi*x)*(1+1/(12*x)+1/(288*x^2)-139/(51840*x^3)+ ...)
so its logarithm will be:
z = ln(x!) = x*ln(x)-x+ln(2*Pi*x)/2+ln(1+1/(12*x)+1/(288*x^2)-139/(51840*x^3))+ ...
and of course your final result is simply e^z. The resulting, non-iterative HP-71B BASIC program is:
10 B=365 @ N=23
20 DISP EXP(FNF(B)-FNF(B-N)-N*LN(B))
30 END
90 DEF FNF(X)=X*LN(X)-X+LN(2*PI*X)/2+LN(1+1/12/X+1/288/X^2-139/51840/X^3)
where DISP is the statement used to show a result on the display (aka PRINT, VIEW, etc), EXP is the exponential function e^x, LN is the natural logarithm function, PI is 3.14159265359, and FNF is a user-defined function, which can be substituted for a subroutine in other languages, passing it
the appropriate parameter. Upon running, it promptly returns:
>RUN
.492702772644
which agrees well with your .493 and is accurate to nearly 8 decimal places (actually, the exact result is 0.49270276567601459277458277+). If still more accuracy is desired, you can use some additional terms to the Stirling's series, have a look at this link:
Stirling's series numerators
and its twin for the denominators. You can get easily 12+ decimals by using more terms. On the other hand, if you don't need that much accuracy, you can save RAM and execution time by deleting terms from the Stirling's series above, like this:
90 DEF FNF(X)=X*LN(X)-X+LN(2*PI*X)/2+LN(1+1/12/X)
though you should test if the accuracy you then get is adequate for your intended purposes, particularly for large N. For your specific numeric example, you still get 8 decimals !
Hope that helps. Best regards from V.
Edited: 17 June 2004, 10:26 a.m.
▼
Posts: 33
Threads: 12
Joined: Jan 1970
This is rich stuff! We're definitely on the right track. I've got a monthly newsletter to get out after which I'll try to digest this. Thank you, Valentin. You just may have solved my problem. I never in my life have heard of a Pochhammer Symbol. It sounds like a novel by Frederick Forsyth.
Posts: 1,792
Threads: 62
Joined: Jan 2005
This problem has an application in the well-known statistical exercise called the "birthday problem".
"Among 25 randomly-selected people in a room (with no twins), what is the probability that at least two people in the room have the same birthday?"
To answer this question, determine the probability of no one sharing a birthday. (Ignore leap year):
prob = (364/365)*(363/365)*...*(341/365)
= Perm(364,24) / 365^24
= 0.43130
So, the chance that two of 25 random people have the same birthday is (1-0.43130) * 100% = 56.87%
Rather counter-intuitive, no?
Of course, for larger values in the quotient, a term-by-term looped approach must be programmed.
-- Karl S.
▼
Posts: 58
Threads: 15
Joined: Mar 2006
When I saw the original posting, it immediately reminded me of the 'birthday problem'. About 20 years ago, I tried it out at our office. We were 20 persons at that time. And yes, there were two with the same birthday. Happened to be a girl who had recently started to work with us and myself.
The counter-intuitive side of it appears because one may tend to think of the chances that two particular persons share birthdays - which is something quite different!
▼
Posts: 75
Threads: 5
Joined: Jan 1970
My first sweetheart's dad pulled that trick at a dinner party when I was in my mid teens (he was a statistician who worked for ICI). The occasion is etched on my mind to this day. I was convinced that professionals that played with numbers were either demigods or magicians--a feeling reinforced when I read some of the postings on this forum.
Suffice it to say that I never could cut the mustard in the numbers game whic is probably why I remain in awe of those that do. ;-)
Cameron
xyzzy
Posts: 33
Threads: 12
Joined: Jan 1970
It is indeed. I've written a very provocative manuscript I'm trying market. I can't tell you all much about it now for the sake of confidentiality but when it's published, everyone is welcome to it.
Posts: 33
Threads: 12
Joined: Jan 1970
Karl, this is the best approach for solving this problem. Not only is it easy to program for a calculator with a SOLVER, it is also easy to explain to groups as to what is going on mathematically.
I'm grateful to you and everyone for your interest in this topic. It has fascinated me for years. I've had one paper published on it and I have a newspaper interested in an update.
I think this is really a great forum. Extremely helpful!
Posts: 2,448
Threads: 90
Joined: Jul 2005
Hi,
I couldn't resist, and so I wrote a loop lat night, both on the 32s and the 11c. The 11c is much slower. The 32s prompts for the number of terms (i) and the starting numerator (A). The 11c must be hand-laoded to I and 0 respectively before GSB A to start. 32s starts with xeq M.
Here is the 11c: hERE IS 32S:
LBL B LBL B
365 365
/ /
RTN RTN
LBL E LBL E
RCL 2 RCL B
RTN RTN
LBL A LBL M
RCL 0 INPUT A
GSB B INPUT i
STO 2 RCL A
1 XEQ B
STO-0 STO B
RCL I 1
X<>Y STO-A
- STO-i
STO I
GTO C GTO C
LBL C LBL C
RCL I RCL i
X=0? X=0?
GTO E GTO E
RCL 0 RCL A
GSB B XEQ B
STO * 2 STO * B
1 1
STO - 0 STO - A
RCL I STO -i
X<>Y GTO C
-
STO I
GTO C
regards,
Bill
Edited: 17 June 2004, 10:29 a.m.
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
Bill --
I just couldn't resist writing a tighter program for the HP-11C/15C/34C for the "birthday problem" that is shorter and simpler in flow, requires no initialization of memory, uses fewer labels, and utilizes the built-in loop-counter function "DSE":
LBL "A"
STO 0 store number of people
X<->Y
INT
STO 1 store integer number of days
LST x
FRAC
STO 2 store fractional number of days
RCL 1
RCL 0
-
EEX
3
/
RCL 1
+
STO I loop counter variable is set
RCL 2
STO + 1 replace days in year with full floating value
1
STO 3 initialize answer variable
LBL "0" start of processing loop
DSE I
GTO "1"
RCL 3
RTN display final result and stop
LBL "1" calculate term
RCL I
INT
RCL 2
+
RCL 1
/
STO x 3
GTO "0"
Usage:
(number of days in year -- may be non-integer)
ENTER
(number of people)
GSB "A"
:
:
(read probability of all unique birthdays)
This program can accommodate leap years by representing 365.25 days (or a more precise value) in a year.
This is somewhat useful, because you can't do
Perm(365.25,25) / (365.25^25)
Edited: 22 June 2004, 2:08 a.m. after one or more responses were posted
▼
Posts: 2,448
Threads: 90
Joined: Jul 2005
Hi Karl,
Thank you very much for your post. I was really hoping to get a critique / response!
I figured I would get a response regarding the possibility of using DSE. The thing was, I can never remember the format for DSE without checking the manual--so I made my own counter on the spot.
I discovered that you can run my 32s version on the 20s---with no modification except the Loop being " / 365 =" instead of "365 /" (There is no DSE/ISG on the 20s).
The other conundrum (my weakest part with loops) was getting the count right---you know--the "when do you count--before or after the action" and so when I tried later to convert to DSE and use a single loop, I ended up with a program that ran fine, except that the counter was off by one.
So thank you for the chance to see how to make it work correctly.
Logic and problem solving are interesting. My solution I posted was very much the first attack---and I programmed it in that classically useful way of doing it manually and recording what you do. It is for that reason that I ended up with an initial solve, store to a result variable, then go into a loop and add to the result variable. I knew it should be possible to do it in one fell swoop---but I couldn't see it right then and there.
But programming is a skill, and you only get better with practice. and to that end, this was a very good exercise.
The problem for me is that I am not a programmer---yet every year, or sometimes after 2 or three years, I find myself needing to write programs---in a very specialized languege----and so while rusty, I do need to know what I am doing. (The language is written for developing stability caclulations for ships---and only infrequently do I need to write a custom program to attack an unusual problem).)
So, thanks again!
Best regards,
Bill
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
... Check my previous post (to which you commented) for an edited version that can calculate correct results with non-integer days in a year.
Posts: 33
Threads: 12
Joined: Jan 1970
OK, my newsletter is mailed so I can get back to the ‘birthday paradox.’ Below is my contribution to the recent spate of programs to solve the birthday paradox. My program appears to be beefy at first glance, but it does a lot of things automatically, including prompts. First, in this problem, what you’re really solving for is the following:
1 ~ [(364/365) * (363/365) * 362/365)…]
So that leaves us with three variables: 365, or what I call ‘max’; number of iterations, which I call ‘group’: and the probability at a certain group count. You would never want to solve for ‘max’ but you have to input it.
This program can possibly be done with fewer steps by taking ‘max’ to the group number power. You can try that; it would be interesting, but I don’t think you’re going to cut down very many steps. I’m not a proficient BASIC programmer but I’m sure a shorter program can be written in either in Virtual or old BASIC. There is no pride of authorship here so have at it. We can all learn.
Here’s the way the program looks on the HP 19BII, using the solver:
1 – (Perm ( X:Y) / (365^Y)) = P where X = universe (365), Y= no. in group (23, and P = target or resultant probability of a matching birthday.
Step Action Explanation for action
01 LBL BDP Label BirthDay Paradox
02 CLRG Clear all registers
03 CF01=GRP Reminder of which variable to solve
04 PROMPT
05 1 Explained in text above
06 STO 3
07 MAX? What is size of universerse?
08 PROMPT
09 STO 2
10 PROB? What is target probability, answer will be in group size
11 PROMPT
12 STO 04
13 GRP? What is group size, answer will be in probability
14 PROMPT
15 STO 05
16 LBL 01 LBL 01 begins the routine to solve for group size
17 1
18 STO +01
19 RCL 02
20 RCL 01
21 ---
22 RCL 02
23 /
24 RCL 03
25 *
26 STO 03
27 1
28 ---
29 CHS
30 STO 06
31 FS? 01 Sets flag to branch to probability routine, otherwise continue
32 GTO 02
33 RCL 04 Recalls target probability to compare with value in Step 30
34 X>Y?
35 GOTO 01 If probability doesn't satisfy, return to LBL 01 & iterate
36 STOP
37 VIEW 01 View the group size since the routine has stopped
38 STOP
39 VIEW 06 View exact probability since it won't exactly = target
40 LBL 02 Starts new routine solvling for Probability
41 RCL 05 Steps 39 & 40 compare the group counter and
42 RCL 01 compares in Step 41
43 X=Y? Here can use X =Y because number in group will be an integer
44 GTO 03 If they equal then branch to Step 44. Read probability
45 GTO 01
46 LBL 03
47 VIEW 06
48 END
|