Issues with real precision

[Fernando Gouvea]

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

[dim...@gmail.com]

here is how to do it exactly, and convert the result to a float.

sage: sage: g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
....: sage: f(n,r)=sum(g(n,k,r),k,0,n)
....: sage: f(50,200).unhold().n()
0.398287031966894
....: sage: f(365,2000).unhold().n()
0.216119451633218


There should be less brute-force ways too,
bu this is the 1st that comes to mind.

HTH
Dima

>
> Fernando
>
> --
> =============================================================
> Fernando Q. Gouveahttp://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
>

> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/fd3261d8-3ee5-48a8-ba6e-40530790b92f%40colby.edu.

[Dima Pasechnik]

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

HTH

Dima

[Emmanuel Charpentier]

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,

[Gareth Ma]

Note that you can wrap it in `Decimal` or `Fraction`, which are both builtin Python libraries.

[Emmanuel Charpentier]

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.21611945163321844

A tad more muscular than using Sage…

Thanks a lot ! I learned something useful today…