The calculated result is mathematically correct.
The code does not sample whole circles, and the deviation from mean is almost entirely due to the final fractional circle. The deviation is a mathematical property of the sequence, not an artifact of floating-point arithmetic.
With lambda = 3,005, the summation samples 3,005 points distributed around the circle. These points are repeated as the iterations go on. With nIterations = 1,000,000, the first 2,340 points are sampled ⌈1,000,000 / 3,005⌉ = 333 times, and the remaining 665 points are sampled 332 times. Given this, it is unsurprising the mean of the samples differs from the mean of complete circles. The magnitude of the differences, about 10−3, is proportional for the number of whole circles completed.
We can verify the errors of floating-point arithmetic are slight by setting nIterations to an integer multiple of lambda. Mathematically, this would produce a sequence with mean 5 (from the mean parameter), so any deviation from 5 is due to imperfect arithmetic. With nIterations set to lambda (3,005), I get a final mean of 5.00000000000000710542735760100185871124267578125 (using code from the question reimplemented in C). With nIterations set to 1000•lambda (3,005,000), I get a final mean of 5.0000000000002788880237858393229544162750244140625.
In the latter cases of nIterations, we see tiny differences that are due to errors in floating-point arithmetic. With the original nIterations, the differences is almost entirely due to the incomplete final circle. The difference is mathematically real, not an artifact of floating-point arithmetic.
Note that the code at the Go playground link differs from the code in the question. It prints the final result using an extra mean(gen()) evaluation, so it effectively performs nIterations+1 iterations. To reproduce the whole-circle results, set nIterations to one less than the desired multiple, such as 3,004 or 3,004,999 instead of 3,005 or 3,005,000.
(Bonus note: the 3,005-iteration result is eight ULP from 5, and the 3,005,000-iteration result is 314 ULP from five.)
Finding the desired mean.
Per additional request from OP, here is code (in C) that finds the mean.
#include <math.h>
static double f(void)
{
static const double π = 0x3.243f6a8885a308d313198a2e03707344ap0;
static const double
λ = 3005,
a = 2,
µ = 5,
ω = 2 * π / λ;
static double t = 1234.5678;
return a*cos(ω * t++) + µ;
}
#include <stdio.h>
/* Given a function f that returns a*cos(ω*t)+µ for unknown constants a, ω,
and µ and a t that starts with an unknown value and increments by one in
each call, calculate µ.
We do this by finding extrema in the function to estimate its period then
sampling points equally spaced around the cycle and applying an algebraic
solution from Maple. (Judging from the simplicity of the solution, I
expect it comes from simple trigonometric properties, but I do not know its
derivation.) Maple was asked to solve:
a*cos(w*(t+0))+u = t1,
a*cos(w*(t+1))+u = t2,
a*cos(w*(t+2))+u = t3, and
a*cos(w*(t+3))+u = t4,
for u, a, w, and t. Since w is an unknown constant, the scale of 0, 1, 2,
and 3 are irrelevant to the solution for u, so it will work for any equally
spaced points in the series produced by the function (barring issues such
as all the points having the same phase so they represent only one point in
the cycle instead of four). Since t is unknown, it will work for any
starting point.
Note that if the wavelength, λ (see sample f above), is an even integer, we
can simply average the two extremes to find the mean. However, λ is
unknown, and the solution below is general.
Mathematically, adjacent samples would suffice. However, testing showed
artifacts due to floating-point precision. So spacing that is a
substantial portion of the period is used. I have not analyzed the
numerical behavior in this situation.
*/
int main(void)
{
// Look for some change in the function.
double t0 = f(), t1;
for (int i = 0; i < 100; ++i)
{
t1 = f();
if (t0 != t1)
break;
}
if (t0 == t1)
printf("Function appears to be constant, mean is %.99g\n", t0);
else
{
// Find first change of direction and then distance to second change.
int count = 0;
if (t0 < t1)
{
do { t0 = t1; t1 = f(); } while (t0 <= t1);
do { ++count; t0 = t1; t1 = f(); } while (t0 >= t1);
}
else
{
do { t0 = t1; t1 = f(); } while (t0 >= t1);
do { ++count; t0 = t1; t1 = f(); } while (t0 <= t1);
}
// Sample function at roughly quarter-circle spacings.
count /= 2;
if (count < 1) count = 1;
for (int i = 1; i < count; ++i) f();
double t2 = f();
for (int i = 1; i < count; ++i) f();
double t3 = f();
for (int i = 1; i < count; ++i) f();
double t4 = f();
// Apply solution from Maple.
printf("u = %.99g.\n",
(t3*t3 + t3*t1 - t2*t2 - t2*t4) / (t1 + 3*t3 - 3*t2 - t4));
}
}
Sample output:
u = 4.99999999999999733546474089962430298328399658203125.
