[
https://issues.apache.org/jira/browse/STATISTICS-100?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=18114044#comment-18114044
]
Alex Herbert commented on STATISTICS-100:
-----------------------------------------
Evaluation of the zeta function can be performed using the Euler-Maclaurin
formula to approximate the summation using a continuous integral. Formulas are
provided in:
{noformat}
Johansson (2015)
Rigorous high-precision computation of the Hurwitz zeta function and its
derivatives
Numerical Algorithms (69) 253–270
{noformat}
The summation is split into two summations, and the second is computed using
the Euler-Maclaurin formula with an integral I, tail T and residual error term
R:
{noformat}
N-1 oo
zeta(s, a) = sum f(k) + sum f(k) = S + I + T + R
k=0 k=N
1
f(k) = ------
s
(a+k)
1-s
,-oo 1 (a+N)
I = | ------ dt = --------
-' N s s-1
(a+t)
/ B (s) \
1 | 1 M 2k 2k-1 |
T = ------ | - + sum ----- --------- |
s | 2 k=1 (2k)! 2k-1 |
(a+N) \ (a+N) /
R = (not evaluated)
B = Bernoulli number
2k
___n-1
(s) = | | (x+i) (rising factorial Pochhammer function)
n | |i=0
{noformat}
Johansson extends these formulas to the complex argument s with real(s) > 1. N
and M must be carefully chosen to achieve the desired precision P. The paper
suggests N ~ M ~ P for P bits of precision. In practice the N and M can be
lower than 53 when evaluating for real s in double precision. I have tested
several implementations with varying s > 1 and a in the range of a positive
integer (required for the domain of the discrete Zipf distribution). Results
will be added below.
> Zipf distribution can use the Hurwitz zeta function to avoid summation over
> all PDF values
> ------------------------------------------------------------------------------------------
>
> Key: STATISTICS-100
> URL: https://issues.apache.org/jira/browse/STATISTICS-100
> Project: Commons Statistics
> Issue Type: Improvement
> Components: distribution
> Affects Versions: 1.3
> Reporter: Alex Herbert
> Assignee: Alex Herbert
> Priority: Minor
>
> The [Zipf distribution
> (Wikipedia)|https://en.wikipedia.org/wiki/Zipf%27s_law#Formal_definition]
> uses a probability mass function (PMF) of:
> {noformat}
> 1 1
> pmf(k; N, s) = -- . ----
> s H
> k N,s
> {noformat}
> The normalising constant H_N,s is the N-th generalised harmonic number of
> order N of s.
> The cumulative probability for range [a, b] is computed by summation of the
> power term k^-s for k in [a, b]. This can be millions of power terms for one
> cumulative probability when the support [1, N] is large.
> When s is above 1 the cumulative probability can use the [Hurwitz zeta
> function (Wikipedia)|https://en.wikipedia.org/wiki/Hurwitz_zeta_function]:
> {noformat}
> oo 1
> zeta(s, a) = sum ------
> k=0 s
> (k+a)
> cdf(k; N, s) = (zeta(s, 1) - zeta(s, k+1)) / H_N,s
> sf(k; N, s) = (zeta(s, k+1) - zeta(s, N+1)) / H_N,s
> H_N,s = zeta(s, 1) - zeta(s, N+1)
> {noformat}
> The Hurwitz zeta function for s > 1 is absolutely convergent. This function
> can be used to significantly improve the performance of the Zipf distribution.
--
This message was sent by Atlassian Jira
(v8.20.10#820010)
