rand_distr: Add Zipf distribution

[vks]

Fixes #1069.

[vks]

cc @kaimast

[vks]

I did not compare to the reference implementation in numpy yet.

[dhardy]

May be worth comparing this? https://jasoncrease.medium.com/rejection-sampling-the-zipf-distribution-6b359792cffa

[vks]

May be worth comparing this? https://jasoncrease.medium.com/rejection-sampling-the-zipf-distribution-6b359792cffa

On a first look, I'm not sure this method is different, but I would have to do the math to be sure.

[vks]

@dhardy The parametrization of the Zipf distribution you linked is different:

• It uses the rank N and the exponent s >= 0, returning an integer from [1, N]. This follows Wikipedia.
• The current implementation only has the exponent a >= 0 as a parameter and returns integers >= 1. This follows numpy.

Which paramtrization should I implement?

[dhardy]

Sorry, I don't feel qualified to answer that.

[vks]

@dhardy Fair enough, I'm also not sure what's best without having a use case in mind.

@kaimast Any preferences? What's your use case for sampling the Zipf distribution?

[vks]

I also found an R library using the N-s parametrization.

The a parametrization is a limit of the N-s parametrization for N -> oo. So maybe N-s is a better choice? Or both should be supported?

Wikipedia calls a the zeta distribution, so maybe we should do the same.

[vks]

@dhardy I chose to implement both distributions. Zeta is equivalent to numpy.random.zipf, Zipf is equivalent to the Java implementation in the blog post you linked. I had to derive the correct equations for the s = 1 special case.

[saona-raimundo]

Hi @vks!
@dhardy asked me if I could help reviewing this implementation, so I hope you don't mind some comments here and there :)

Looking at the previous comments you had quite a journey with names (and yeah, it is not very clear). But you got it right!

I have looked only at the implementation of Zeta for now. I added some suggestions on the documentation, but wanted to highlight that I have some questions about the sampling algorithm: I found the reference of the algorithm [Devroye86], so we have proved correctedness :) and my only doubts is about floating point representability... Well, we will talk more in the comments, I think.

Naming

I had a look at some references and leave a summary here. Maybe this is useful for the documentation, but I will leave it to you.

1. The general family of distributions is named Lerch [JKK05], after the Lerch zeta function. It is also called Hurwitz-Lerch [USB18].

2. There are many related distributions (at least 8 variants) [JKK05].

3. As you noticed, there are two variants (sometimes called the same): Zeta and Zipf.

• Zeta distribution has support over all integers 1, 2, ...
• Zipf distribution has support only over 1, 2, ..., N (it is simply the truncated version).

4. Other names for Zeta are: discrete Pareto, Riemann-Zeta, Zipf and Zipf–Estoup.

5. Special cases of Zipf distribution have names too: Estoup for s = 1 and Lokta for s = 2.

References

• [JKK05] Section 11.2.20, page 526.
• [Devroye86] Section 6.1, page 551.
• [USB18] Chapter 3, Section 3.2.3, page 117.

[vks]

I would suggest saying something in the documentation about returning floats for integer random variables, but that deserves its own discussion and PR. This would involve investigating what is the effective domain of our implementations... information which is valuable for the user, but hard to get (or depends on parameters).

There was some discussion in #987 and #1093, which resulted in some documentation in the Rand book.

[saona-raimundo]

These are the comments about Zipf.

I have some disagreements with some of the computations, I tried to explain the reasoning in each case. Feel free to ask for more clarification.

Notably, I also propose to introduce one more error variant NTooBig.

[saona-raimundo]

There was some discussion in #987 and #1093, which resulted in some documentation in the Rand book.

I see! Thanks, I did not know about it :)

[vks]

Another case we might want to consider: For large a, the parameter b will be infinite, which results in all x being accepted. This is probably the best we can do, but we could also rather return an error when Zeta is constructed.

[vks]

Because this PR is already quite large, I would like to leave the extended distribution tests for a follow-up PR.

[saona-raimundo]

@dhardy we are finished with the review and agree on adding distributional tests in a subsequent PR.

[dhardy]

Great. On that basis I approve merging. And thanks for stepping in to review.

[vks]

@saona-raimundo Thanks for the extensive review!

@dhardy Should I squash before merging?

[dhardy]

@vks can't say I care too much whether or not it gets squashed. The reason I didn't merge myself is because I wasn't quite sure whether you were ready to.