85 views

Skip to first unread message

Feb 23, 2024, 5:00:50 PMFeb 23

to sage-support

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

Feb 23, 2024, 5:16:37 PMFeb 23

to sage-s...@googlegroups.com

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

>

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

Feb 23, 2024, 5:23:20 PMFeb 23

to sage-support

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

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

Feb 24, 2024, 8:54:51 AMFeb 24

to sage-support

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,

Feb 24, 2024, 12:11:25 PMFeb 24

to sage-support

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

Feb 24, 2024, 5:36:34 PMFeb 24

to sage-support

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) :

A tad more muscular than using Sage…

Thanks a lot ! I learned something useful today…

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu