Cut the Cards

07302020, 08:00 PM
Post: #1




Cut the Cards
PPC Journal V5N4 Page 13 poses a question about a deck of cards. Choose a card at random from a deck of 52 and put the card back. What is the average number of choices required to see all 52 cards?
I simulated the experiment in C++ and got the right answer (236) but I can't figure how to do this analytically. Is it possible? If you simplify the problem down to a deck of 2 cards then it's easier: you choose the first card, and then it's just a matter of the average number of choices to get the second: 1*1/2 + 2*1/4 + 3*1/8 ..., but can that be extended to more cards? Dave 

07302020, 08:58 PM
Post: #2




RE: Cut the Cards
(07302020 08:00 PM)David Hayden Wrote: Choose a card at random from a deck of 52 and put the card back. First card picked must never been seen, thus only 1 try needed to get 1st card. Probability of getting the next "new" card = 51/52, or averaged 52/51 tries to get 2nd card ... We can get averaged picks either way: XCas> sum(1/t, t=1 .. 52) * 52.0 → 235.978285436 XCas> int((t^521)/(t1), t=0 .. 1) * 52.0 → 235.978285436 → averaged 236 picks to see all 52 cards 

07302020, 09:49 PM
Post: #3




RE: Cut the Cards
So, it approaches x*(ln(x)+0.57721), where the latter constant is the EulerMascheroni constant. This isn't as accurate as doing it analytically but is good for fast calculations. For instance, if you picked a random second of a year every second, how many years would it take? About 17.84 total.
So many signals, so little bandwidth! 

07302020, 10:21 PM
Post: #4




RE: Cut the Cards
Here is simulations, using Lua
Note: math.random(n) return 1 .. n, each with same probability. Code: function run(tests,cards)  averaged picks lua> run(1e6, 52) 235.970842 

07312020, 12:24 AM
(This post was last modified: 07312020 12:51 AM by John Keith.)
Post: #5




RE: Cut the Cards
(07302020 09:49 PM)Jim Horn Wrote: So, it approaches x*(ln(x)+0.57721), where the latter constant is the EulerMascheroni constant. Or in other words, n*H(n), where H(n) is the nth harmonic number. (Albert's first formula in post #2) (07302020 09:49 PM)Jim Horn Wrote: This isn't as accurate as doing it analytically but is good for fast calculations. Actually it is very accurate as the first image in the Wikipedia link shows: LN(n)+gamma is the asymptotic limit of H(n) as n approaches infinity. 

08212020, 11:00 PM
Post: #6




RE: Cut the Cards
(07302020 08:58 PM)Albert Chan Wrote: First card picked must never been seen, thus only 1 try needed to get 1st card. I might jumped too quickly to say averaged 1/p tries to get the next card. Let probability of getting the card = p ; of not getting it = q = 1p Expected tries = p + 2qp + 3q²p + 4q³p + ... = (1q) (1 + 2q + 3q² + 4q³ + ...) = (1 + 2q + 3q² + 4q³ + ...)  (q + 2q² + 3q³ + ...) = 1 + q + q² + q³ + ... = 1/(1q) = 1/p For n > 0, we have: \(\Large {1 \over (1x)^n}\normalsize = \displaystyle \sum_{k=0}^∞ \binom{n}{k}(x)^k = \sum_{k=0}^∞ \binom{k+n1}{n1} x^k \) Examples: XCas> series(1/(1x)^2, x, polynom) → 1+2*x + 3*x^2 + 4*x^3 + 5*x^4 + 6*x^5 XCas> series(1/(1x)^3, x, polynom) → 1+3*x + 6*x^2+10*x^3+15*x^4+21*x^5 // triangular numbers coefficients XCas> series(1/(1x)^4, x, polynom) → 1+4*x+10*x^2+20*x^3+35*x^4+56*x^5 // tetrahedral numbers coefficients 

08242020, 01:57 PM
(This post was last modified: 08242020 02:07 PM by Gerson W. Barbosa.)
Post: #7




RE: Cut the Cards
(07312020 12:24 AM)John Keith Wrote:(07302020 09:49 PM)Jim Horn Wrote: So, it approaches x*(ln(x)+0.57721), where the latter constant is the EulerMascheroni constant. The H(n) function can be easily defined on the HP 50g: « 1. + Psi 1. Psi  » 'H' STO For instance, 52 ENTER ENTER H × → 235.978285436 

08242020, 09:49 PM
Post: #8




RE: Cut the Cards
Thanks for all those formulas and explanations.
Here is an implementation of the brute force algorithm for the experimental analysis: Code:
Usage: CARDS(numberofcards, numberofanalysis) Example: CARDS(52, 1000) Returns: 235.951 Thibault  not collector but in love with the few HP models I own  Also musician : http://walruspark.co 

08252020, 06:14 PM
(This post was last modified: 08282020 01:23 AM by Albert Chan.)
Post: #9




RE: Cut the Cards
(07302020 09:49 PM)Jim Horn Wrote: So, it approaches x*(ln(x)+0.57721), where the latter constant is the EulerMascheroni constant. Replacing ln(x) by ln(x+0.5) gives much better estimate for ψ(x+1). ψ(x+1) = slope of ln(gamma(t)) at t = x+1. When x is large: \(\psi(x+1) ≈ {\ln\Gamma(x+2)\;\;\ln\Gamma(x)\over 2}=\ln\sqrt{(x+1)(x)} ≈ \ln(x+0.5)\) >>> x = 52 >>> H = log(x+0.5) + 0.57722 >>> print format(x*H, 'g') 235.978 see https://en.wikipedia.org/wiki/Digamma_fu...roximation 

08252020, 11:41 PM
(This post was last modified: 08262020 12:00 AM by Gerson W. Barbosa.)
Post: #10




RE: Cut the Cards
FWIW,
« 1. Psi + EXP 6. OVER INV 12. + / + 1.  » 'H¹' STO Example: 1100 H → 7.58073560027 H¹ → 1100 H¹ should be accurate to 11 or 12 digits for arguments equal or greater than 7.5 

08262020, 03:06 AM
(This post was last modified: 09012020 02:22 PM by Albert Chan.)
Post: #11




RE: Cut the Cards
(08252020 11:41 PM)Gerson W. Barbosa Wrote: FWIW, Here is another version, using asymptotic estimate: \(\quad\exp(\psi(x+0.5)) ≈ x + \large{1\over24x} \) Let k = exp(ψ(x+0.5)), n = x  0.5 → k = exp(H(n) + ψ(1)) ≈ x + 1/(24x) If x is big, x ≈ k, we have: \(n ≈ k  \large{1\over 24k}  {1\over2}\) HP50g inverse harmonic function: « 1. Psi + EXP DUP 24. * INV  .5  » Both formulas perform well for large H, but this new version is also good for smaller H >>> from mpmath import exp, euler, harmonic >>> f1 = lambda k: k + 6/(1/k+12)  1 >>> f2 = lambda k: k  1/(24*k)  .5 # new version >>> for n in range(50,100,10): ... h = harmonic(n) ... k = exp(h  euler) ... e1, e2 = n  f1(k), n  f2(k) # underestimated errors ... print '%d\t%.10f\t%.10f\t%g' % (n, e1, e2, abs(e1/e2)) ... 50 0.0000013228 0.0000000364 36.3546 60 0.0000009261 0.0000000212 43.7606 70 0.0000006844 0.0000000134 51.1669 80 0.0000005263 0.0000000090 58.5736 90 0.0000004172 0.0000000063 65.9804 

08262020, 08:23 AM
Post: #12




RE: Cut the Cards
(08262020 03:06 AM)Albert Chan 0.247866177136� 54 Wrote: HP50g inverse harmonic function: « 1. Psi + EXP DUP 24. * INV  .5  » For H less than 2.0701076 the first formula is slightly better: « DUPDUP H H¹  SWAP DUP H H¹2  / ABS » 50 → 36.4435261708 10 → 6.75851700917 2.0701076 → 1. 1 → 0.290013346857 0.5131 → 2.30966647121E5 0.5 → 7.2984301462E3 0.2 → 0.162258684815 0 → 0.247866177136 

08262020, 02:13 PM
(This post was last modified: 08312020 08:31 PM by Albert Chan.)
Post: #13




RE: Cut the Cards
I added a tiny correction to compensate errors for smaller H.
HP50g H^{1} : « 1. Psi + EXP DUPDUP SQ 24. * 2.7 + /  .5  » Code: from mpmath import exp, euler, harmonic, psi >>> test(range(10) + range(10,60,10)) n err(f1) err(f2) err(f3) 0 0.0031607566 0.0127518672 0.0067666262 1 0.0003177730 0.0010957186 0.0001620997 2 0.0002588484 0.0002719596 0.0000171164 3 0.0001712995 0.0001037256 0.0000035307 4 0.0001178786 0.0000497922 0.0000010529 5 0.0000852046 0.0000275594 0.0000003957 6 0.0000641850 0.0000167992 0.0000001742 7 0.0000499820 0.0000109785 0.0000000859 8 0.0000399758 0.0000075616 0.0000000463 9 0.0000326776 0.0000054263 0.0000000266 10 0.0000271983 0.0000040243 0.0000000162 20 0.0000076840 0.0000005432 0.0000000006 30 0.0000035570 0.0000001651 0.0000000001 40 0.0000020419 0.0000000705 0.0000000000 50 0.0000013228 0.0000000364 0.0000000000 If needed, we could apply Newton's method to get better estimate. >>> h = harmonic(10) # h = 2.9289682539682538 >>> x = f3(exp(heuler)) # x = 10.000000016205583 >>> x  (harmonic(x)  h) / psi(1, x+1) mpf('10.0') 

08262020, 06:13 PM
Post: #14




RE: Cut the Cards
D’oh, my latest improvement ended up being equivalent to your f2 program, same number of bytes as yours (56.5).
That’s an old attempt at exact results, but I wasn’t pleased with it. Perhaps I can make it shorter some day: https://hpmuseum.org/forum/thread6157post55046.html 

08272020, 10:07 PM
Post: #15




RE: Cut the Cards
How about « 1. Psi + EXP DUP 1. + 4. * INV OVER 1. EXP / + INV OVER 24. * + INV  .5  » ?
Again, best for smaller H. A previous attempt ( « 1. Psi + EXP DUP SQ 24. * 1. EXP + INV NEG 1. + * .5  » ) gave the same results as yours when I replaced 1 EXP ( e ) with 2.7. 

08282020, 09:26 PM
(This post was last modified: 08292020 02:51 AM by Albert Chan.)
Post: #16




RE: Cut the Cards
(08262020 02:13 PM)Albert Chan Wrote: I added a tiny correction to compensate errors for smaller H. (08272020 10:07 PM)Gerson W. Barbosa Wrote: How about « 1. Psi + EXP DUP 1. + 4. * INV OVER 1. EXP / + INV OVER 24. * + INV  .5  » ? Just to be clear, both corrections are not based on actual continued fraction of e^ψ(x+1/2) My correction (and probably Gerson's version) are based on curve fitting, for small H. I picked 2.7 instead of e to get better accuracy at bigger H, it is simply a tradeoff. To build e^ψ(x+1/2) continued fraction with XCas: XCas> k := exp(Psi(x+0.5)) XCas> limit(k/x, x=inf) → 1 XCas> limit(1/(k x)/x, x=inf) → 24 XCas> limit(1/(1/(k x) 24x)/x , x=inf) → 10/37 XCas> limit(1/(1/(1/(k x) 24x) 10/37*x)/x , x=inf) → 0 ??? The last one is likely a XCas bug, Mathematica get the coefficient: 689976/74381 \(\qquad\qquad\exp( \psi(x+1/2)) = x +\frac{1}{24 \cdot x +\Large\frac{1}{\frac{10}{37} \cdot x +\frac{1}{\frac{689976}{74381} \cdot x \; +\; \cdots}}} \) To confirm, lets try H^{1}(H(20)): XCas> f(x) := 1/(24*x + 1/(10/37*x + 1/(689976/74381*x))) XCas> h := float(harmonic(20)) → 3.59773965714 XCas> k := exp(h  euler_gamma) → 20.5020317757 XCas> fsolve(x+f(x) = k, x=k)  0.5 → 20.0 If f(x) is small, x ≈ k. We can also iterate for x, like this: XCas> x0 := k XCas> x0 := k  f(x0) → 20.5000002012 XCas> x0 := k  f(x0) → 20.5 XCas> n := x0  0.5 → 20.0 

08282020, 11:39 PM
Post: #17




RE: Cut the Cards
For the inverse digamma function I’ve been using a continued fraction with increasing number of terms. I am sure about only the first half of them, though. It’s difficult to be sure about the other half when one is limited to only 12 digits.
\(\rm{e}^{{x}}+\frac{1}{2+\frac{1}{6\rm{e}^{{x}}\frac{1}{2+\frac{1}{\rm{e}^{{x}1}\frac{1}{2+\frac{1}{\rm{e}^{{x}}+\frac{1}{2}}}}}}}\) 

08292020, 04:02 PM
Post: #18




RE: Cut the Cards
(08272020 10:07 PM)Gerson W. Barbosa Wrote: How about « 1. Psi + EXP DUP 1. + 4. * INV OVER 1. EXP / + INV OVER 24. * + INV  .5  » ? FYI, above H^{1} is equvialent to Ψ^{1} (previous post) Since H(n) = Ψ(n+1)  Ψ(1), we flip the sign of 0.5 to +0.5 XCas> invPsi1(k) := k  1/(24k + 1/(k/e + 1/(4k+4))) + 1/2 \(k \rightarrow k\frac{1}{24\cdot k+ \large\frac{1}{\frac{k}{e}+\frac{1}{4\cdot k+4}}}+\frac{1}{2}\) XCas> invPsi2(k) := k + 1/(2 + 1/(6k  1/(2 + 1/(k/e  1/(2 + 1/(k+1/2)))))) \(k \rightarrow k+\frac{1}{2+\frac{1}{6\cdot k\frac{1}{2+\frac{1}{\frac{k}{e}\frac{1}{2+\frac{1}{k+\frac{1}{2}}}}}}}\) XCas> simplify(invPsi1(k)  invPsi2(k)) → 0 Let's try a few x, to see if it roundtrip. XCas> invPsi(x) := invPsi1(exp^float(x)) // Note: k = e^x XCas> map(range(1,6), x>[x, Psi(invPsi(x))]) \(\qquad\left(\begin{array}{cc} 1 & 0.999999666188 \\ 2 & 2.00000000607 \\ 3 & 3.00000000018 \\ 4 & 4.0 \\ 5 & 5.0 \end{array}\right)\) 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)