In an introductory probability class, one computes the probability of getting all of n possible coupons in r individual purchases. The naive approach with inclusion-exclusion leads to the awful formula
f(n,r) = \sum_{k=0}^n \binom{n,k} \frac{(n-k)^r}{n^r}
Computing this in GP is straightforward:
(10:32) gp > g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r %3 = (n,k,r)->(-1)^(k)*binomial(n,k)*(n-k)^r/n^r (16:43) gp > f(n,r)=sum(k=0,n,1.0*g(n,k,r)) %4 = (n,r)->sum(k=0,n,1.0*g(n,k,r)) (16:43) gp > f(6,12) %5 = 0.43781568062117902081322291656082236787 (16:43) gp > f(50,100) %6 = 0.00016616318861823418331310572392871711604 (16:45) gp > f(50,200) %7 = 0.39828703196689437232172533821416865689 (16:47) gp > f(365,2236)
Doing it in SageMath seems to be much more delicate. This, for example, fails:
sage: g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
sage: f(n,r)=sum(1.0*g(n,k,r),k,0,n)
sage: f(365,2000).unhold()
0.0
It works better if I rewrite the function with ((n-k)/n)^r, of
course. I presume the issue is that the default precision for real
numbers is not enough to handle the computation. How do I tell
SageMath to use a higher precision?
Fernando
-- ============================================================= Fernando Q. Gouvea http://www.colby.edu/~fqgouvea Carter Professor of Mathematics Dept. of Mathematics Colby College 5836 Mayflower Hill Waterville, ME 04901 Reached by phone, a Westbury woman who identified herself as Onorato's granddaughter said he was a retired businessman. He lived alone and never married or had children, she said. -- New York Newsday, December 27, 2003
Le vendredi 23 février 2024 à 23:23:20 UTC+1, Dima Pasechnik a écrit :
[ Snip…]
the normal Python way, without any symbolic sum, would be like this:sage: sage: g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
....: sage: def f(n,r): return math.fsum([1.0*g(n,k,r) for k in range(n+1)])
....: sage: f(365,2000)
0.21611945163321847
This works in Sagemath. It wouldn’t work in Python : the range of the magnitudes of the terms of the alternating sum are way too large for the precision of Python’s floats, necessary if you want to use math.comb. Programming this in Python would need some serious analytical work, or using a multiple-precision integer library, which Sage does for you…
[ Re-snip... ]
HTH,
Le samedi 24 février 2024 à 18:11:25 UTC+1, Gareth Ma a écrit :
Note that you can wrap it in `Decimal` or `Fraction`, which are both builtin Python libraries.
I stand corrected. Python ints are bignums, not 32-bits integers. and fractions.Fractions are rationals. Ordinary integer arithmetic wroks with these rationals. We just need a binomial coefficient function returning an int (math.comb returns a float) :
>>> from math import factorial >>> from fractions import Fraction >>> def ibin(n, k): return Fraction(factorial(n), factorial(k)*factorial(n-k)) if n>=k else 0 ... >>> def g(n, k, r): return (-1)**k*ibin(n,k)*Fraction((n-k)**r, n**r) ... >>> def f(n, r): return float(sum([g(n,k,r) for k in range(n+1)])) ... >>> f(365, 2000) 0.21611945163321844A tad more muscular than using Sage…
Thanks a lot ! I learned something useful today…