[
https://issues.apache.org/jira/browse/STATISTICS-100?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=18114130#comment-18114130
]
Alex Herbert commented on STATISTICS-100:
-----------------------------------------
There are some variations on how to implement the sums.
The sum of the power series in S can be ascending k in [0, N) or descending k
in (N, 0]. The terms reduce in magnitude when k is larger and so the descending
k series should be more accurate by summing small to large terms.
The power term in T can use Math.pow for the power function
{{{}(a+N){^}(2k+1){^}{}}}. Since (a+N) is constant It can use the recurrence
relation: {{{}(a+N)-(x+1) = (a+N)^-x / (a+N){}}}.
The terms can be added in extended precision (e.g. using a double-double sum).
I tested the following implementations:
||Method||Notes||
|zeta1|Ascending k sum; Math.pow for the tail series|
|zeta2|Descending k sum; Math.pow for the tail series|
|zeta3|Descending k sum; recycles (a+N)^-s in the tail series|
|zeta4|Ascending k sum; recycles (a+N)^-s in the tail series; all summation in
extended precision|
I tested using 50 s values in [1, 25) with at 8 values less than 1.01; and 50 a
values in [1, 2^31) with 10 values in [1, 5). Total references values 2500. The
reference used Matlab symbolic toolbox hurwitzZeta function to 30 digits. Data
can also be generated using Python's mpmath zeta function. These agree to at
least 25 digits which is suitable to cross validate the result in double
precision.
The results vary by JVM as the implementation of the Math.pow function is JVM
specific. Here are typical results from JDK Temurin 25.472-b08.
||Method||N||Max||RMS||
|zeta1|9|3.739011012156354|0.6548846555281865|
|zeta2|9|3.698235736226918|0.7690632807907499|
|zeta3|5|19.5732122922785|2.1385715706164765|
|zeta3|6|3.739011012156354|0.6039110504123638|
|zeta3|7|3.698235736226918|0.5979147725867169|
|zeta3|8|3.463160912553975|0.5808808049509839|
|zeta3|9|3.739011012156354|0.6008010722987028|
|zeta3|10|3.463160912553975|0.5854149351464514|
|zeta3|11|3.78356780754505|0.5927549049183026|
|zeta3|12|3.78356780754505|0.5921653127320653|
|zeta4|9|3.463160912553975|0.49242628834265467|
zeta3 has a lower RMS than the other double precision implementations. zeta4 is
best for RMS error. The max error is unchanged and may be limited by the double
precision Math.pow function. zeta4 RMS is only lower than zeta3 if the full
double-double precision is used throughout. Skipping it in the tail series
increases the RMS.
The choice of N takes some tuning. If too low then precision is impacted. Too
high then the sum includes more terms, each of which can have 1 ULP error from
the Math.pow function and this can accumulate into the result.
For the zeta3 N=8 impementation the tail series converges after various M. Here
is a table of M when the tail series is summed to below 1 ULP of the current
sum (S+I):
||Method||N||M||Count||
|zeta3|8|1|385|
|zeta3|8|2|596|
|zeta3|8|3|199|
|zeta3|8|4|129|
|zeta3|8|5|215|
|zeta3|8|6|172|
|zeta3|8|7|223|
|zeta3|8|8|186|
|zeta3|8|9|295|
> 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)
