[
https://issues.apache.org/jira/browse/STATISTICS-101?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=18115618#comment-18115618
]
Alex Herbert commented on STATISTICS-101:
-----------------------------------------
h2. Implementation Details
The zeta distribution is defined over all positive integers. However the
implementation is limited to a 32-bit signed integer by the
DiscreteDistribution interface. This leads to parameterisations that can
truncate the distribution. This occurs as s -> 1.
||s||x||sf(x; s)||
|1.20889000...|2147483647|0.01|
|1.10440113...|2147483647|0.1|
|1.03141836...|2147483647|0.5|
|1.00477515...|2147483647|0.9|
The survival function (sf) can be implemented using the Hurwitz zeta function
for any parameters x and s with efficient evaluation and trivial loss of
precision. But the cumulative probability function (cdf) can be evaluated
either using 1 - sf, or a direct summation of the power terms of the series.
When the sf > 0.5 then the cdf will lose precision. The number of bits of
precision lost is expected to be b-bits for a sf value > 1 - b^(b+1). For
example
||sf||bits lost in 1 - sf||
|0.5|0|
|0.75|1|
|0.875|2|
If the survival function is above 0.5 evaluating the cdf by summation of terms
is increasingly expensive as s -> 1. Allowing 2-bits of precision loss is a
limit of 0.875 with the following number of terms x for different s:
||s||x||sf(x; s)||
|1.04565582...|10|0.875|
|1.02575353...|100|0.875|
|1.01784271...|1000|0.875|
|1.01364468...|10000|0.875|
|1.01104559...|100000|0.875|
|1.00927825...|1000000|0.875|
|1.00799849...|10000000|0.875|
There is always a parameterisation where summation of the terms will be too
expensive for some use cases. Likewise there is no small number of terms where
precision loss can be limited in all cases. The following shows the value of s
for the sf(x=10). Each successive entry in the table will lose another bit of
precision if the cdf = 1 - sf:
||s||x||sf(x; s)||
|1.23838860...|10|0.5|
|1.09852360...|10|0.75|
|1.04565582...|10|0.875|
|1.02204999...|10|0.9375|
|1.01084329...|10|0.96875|
Note that the implementation of the zeta function requires 9 calls to Math.pow
and some additional overhead. So if the number of terms is small it is faster
to compute the sum of the series of power terms rather than the zeta function.
To avoid excessive runtimes I have implemented the distribution using the
survival function as the definitive cumulative probability result. This is
accurate over the entire domain. The cdf is implemented using a sum of the
power series if the number of terms is 10 or less; otherwise it uses 1 - sf.
From the table above this should allow the cdf to maintain accuracy for s >
1.045.
The cumulative probability range function, p(a < x <= b), is implemented in the
same way. This must be considered increasingly less precise as the computed
probability drops below 0.5 and the range b - a >= 10.
h2. High Precision CDF
The cumulative probability can be computed manually using a sum of the power
terms. This must be normalised by the value of the zeta function for s. Note
that the probability for x=1 is the term 1 / zeta(s). So the power terms can be
computed and the distribution used to provide the normalisation factor. Summing
the terms in ascending order of magnitude will reduce cumulative error and can
be done for example using:
{code:java}
double s = ...
ZetaDistribution dist = ZetaDistribution.of(s);
double cdfX = IntStream.range(0, x)
.mapToDouble(y -> Math.pow(x - y, -s))
.sum() * dist.probability(1);
{code}
Summing in descending order of magnitude for the power terms (k=1 to x) can be
used to compute the CDF for all 1 <= k <= x in a single pass loop. There may be
loss of precision for large k due to error in the cumulative sum. This could be
mitigated using an extended precision sum, for example using a double-double
(DD) number from Commons Numbers:
{code:java}
doube s = ...
int x = ...
double[] cdf = new double[x];
double p1 = ZetaDistribution.of(s).probability(1);
cdf[0] = p1;
DD sum = DD.ONE;
for (int k = 2; k <= x; k++) {
sum = sum.add(Math.pow(k, -s));
cdf[k - 1] = sum.doubleValue() * p1;
}
{code}
h2. Sampling
Sampling can be done by inversion of the survival function for probability p.
This will be slow due to the use of a search for x where sf(x; s) - p == 0. A
single call to the survival function is also moderately expensive.
No zeta distribution sampler exists in Commons RNG. But a rejection method is
supplied in:
{noformat}
Devroye, L (2015)
Non-uniform random variate generation.
Springer New York, NY. pp 550-552.
{noformat}
I have implemented this method and tested it where the quantiles of the
distribution are spread over a reasonable range for x. This requires an s
parameter in [1.1, 1.35]. When s -> 1 the samples are very skewed to x=1; when
s -> large the samples are skewed towards infinity which is clipped to
2147483647.
This implementation can be migrated to Commons RNG.
> Add a zeta distribution
> -----------------------
>
> Key: STATISTICS-101
> URL: https://issues.apache.org/jira/browse/STATISTICS-101
> Project: Commons Statistics
> Issue Type: New Feature
> Components: distribution
> Affects Versions: 1.3
> Reporter: Alex Herbert
> Priority: Minor
>
> The [zeta distribution
> (wikipedia)|https://en.wikipedia.org/wiki/Zeta_distribution] is the large N
> limit of the zipf distribution:
> {noformat}
> pmf(x) = k^-s / zeta(s)
> {noformat}
> Where zeta(s) is the Riemann zeta function.
> This can be implemented using the Hurwitz zeta function added for the zipf
> distribution in [STATISTICS-100].
--
This message was sent by Atlassian Jira
(v8.20.10#820010)
