Calculate mean and standard error of a generator in Python

The goal of the following code is to calculate the mean and the standard error of vectors of randomly-generated numbers. I am looking for feedback both on the correctness of the calculation, and on its efficiency.

import numpy as np

def mean_and_stderr(num_of_iterations:int, instance_generator) -> (float,float):
"""
Calculate the mean and standard error of the given generator.

:param instance_generator: a function that accepts no parameters,
and returns either a float or a numpy array.
:param num_of_iterations: number of times to run the instance_generator.

:return a tuple: (mean, standard_error)

Test on a degenerate (constant) generator of numbers:
>>> generator = lambda: 5
>>> mean_and_stderr(100, generator)
(5.0, 0.0)

Test on a degenerate (constant) generator of vectors:
>>> generator = lambda: np.array([1,2,3])
>>> mean_and_stderr(100, generator)
(array([ 1.,  2.,  3.]), array([ 0.,  0.,  0.]))
"""
sum = sumSquares = None
for i in range(num_of_iterations):
x_i = generator()
if sum is None:
sum = x_i
sumSquares = (x_i*x_i)
else:
sum += x_i
sumSquares += (x_i * x_i)
mean = sum / num_of_iterations
variance = sumSquares / num_of_iterations - (mean*mean)
stderr = np.sqrt(variance) / num_of_iterations
return (mean,stderr)

if __name__=="__main__":
generate_uniformly_random_number = np.random.random
print(mean_and_stderr(10, generate_uniformly_random_number))
# Typical output: (0.5863703739913031, 0.026898107452102943)
print(mean_and_stderr(1000, generate_uniformly_random_number))
# Typical output: (0.514204422858358, 0.0002934476865378269)

generate_uniformly_random_vector = lambda: np.random.random(3)
print(mean_and_stderr(10, generate_uniformly_random_vector))
# Typical output: (array([ 0.53731682,  0.6284966 ,  0.48811251]), array([ 0.02897111,  0.0262977 ,  0.03192519]))
print(mean_and_stderr(1000, generate_uniformly_random_vector))
# Typical output: (array([ 0.50520085,  0.49944188,  0.50034895]), array([ 0.00028528,  0.00028707,  0.00029089]))


2 Answers

Formulas are flawed

Both these formulas are incorrect. They seem made up variants of the correct ones.

 variance = sumSquares / num_of_iterations - (mean*mean)
stderr = np.sqrt(variance) / num_of_iterations


Here are the correct formulas:

$$\newcommand{smallm}[0]{\overset{n\ \gg\ m}{\longrightarrow}} \begin{array}{|l|c|c|} \hline \\ & \textrm{Formula} \\ \hline \\ \textrm{Variance} & \sigma² = \frac{1}{N}\sum_{i=1}^N(x_i - \mu)^2 \\ \hline \\ \textrm{Standard Deviation} & \sigma = \sqrt{\sigma²} \\ \hline \\ \textrm{Standard Error} & {\sigma}_\bar{x}\ = \frac{\sigma}{\sqrt{n}} \\ \hline \end{array}$$

Verification

Tool used to verify: Online Standard Error Calculator

Note that the online test is compliant to a sample space, not an entire population. This means the formula used for variance is slightly different to take into account outliers:

$$s² = \frac{1}{n-1}\sum_{i=1}^n(x_i - \bar{x})^2$$

Let's take a fixed input array to test your formulas.

 input array: { 0, 0, 1, 2, 3, 12 }


$$\newcommand{smallm}[0]{\overset{n\ \gg\ m}{\longrightarrow}} \begin{array}{|l|c|c|} \hline \\ & \textrm{Mean} & \textrm{Standard Error} \\ \hline \\ \textrm{OP} & \color{blue}{3} & \color{red}{0.693888666488711} \\ \hline \\ \textrm{Corrected (sample space)} & \color{blue}{3} & \color{blue}{1.69967317119759} \\ \hline \\ \textrm{Corrected (population)} & \color{blue}{3} & \color{blue}{1.8618986725} \\ \hline \end{array}$$

• Thanks a lot! You saved me from a severe bug. I did not use the standard formulas since they require to do two passes on the data: one to calculate the mean $\mu$, and one to calculate the variance $\sigma^2$. Now, I found a way to do the correct calculation in a single pass: mean = sum / num_of_iterations; population_variance = sumSquares / num_of_iterations - (mean * mean); sample_variance = population_variance*num_of_iterations/(num_of_iterations-1); stderr = np.sqrt(sample_variance / num_of_iterations); – Erel Segal-Halevi Jul 13 at 19:01
• Glad to have been of assistance: this MathJax library was not easy to use :) For future questions, I would suggest you add unit tests to your code. This way you can find bugs like this from the get-go. – dfhwze Jul 13 at 19:07
• @ErelSegal-Halevi Btw, read the answers on this question for a single pass calculation of both mean and variance: codereview.stackexchange.com/questions/185450/… – dfhwze Jul 13 at 19:41

I'd want to address the generator, to begin. :param generator: accepts no parameters and returns a vector is false. The generator returns a number (with floating point or not). This is pretty confusing. When I read the code at first, I thought the generator would be an iterable that returns numbers (like range(0,10) for example).

In that case you wouldn't need to pass both the parameters num_of_iterations and generators. Otherwise, I don't think the parameter should be named generator or maybe your documentation should be stronger.

Next thing, you shouldn't initialize sum = sumSquares = None this way. They are numbers, initialize them at zero, that would give you the opportunity to remove your if/else.

sum = sumSquares = 0
for i in range(num_of_iterations):
x_i = generator()
sum += x_i
sumSquares = (x_i*x_i)


Apart from that, the coding style is a little off. Sometimes you use camelCase and sometimes snake_case.

If we were to have an iterator instead of the generator, you could do something like this :

def mean_and_stderr(iterator) -> (float,float):
sum = sumSquares = 0
for x in iterator:
sum = x
sumSquares = (x*x)

mean = sum / num_of_iterations
variance = sumSquares / num_of_iterations - (mean*mean)
stderr = np.sqrt(variance) / num_of_iterations
return (mean,stderr)

• I think whatever_that_is_case is called snake_case (en.wikipedia.org/wiki/Snake_case). – AlexV Jul 12 at 16:33
• @AlexV thanks that's fixed – IEatBagels Jul 12 at 16:34
• It seems like the Python style guide itself refers to that style as lower_case_with_underscores. But it doesn't really matter :-D – AlexV Jul 12 at 16:46
• Thanks for the comments. Probably my use of the word "generator" was confusing. I meant generator in this sense: en.cppreference.com/w/cpp/algorithm/generate i.e., a function without parameters. Like in the phrase "random number generator". Can you suggest a less ambiguous term? Also, the reason I initialize to None is that my "generator" can return not only a single number, but also a numpy array. See the second example in the doctest and the second example in the main program. Is there a more efficient way to handle both numbers and arrays in the same function? – Erel Segal-Halevi Jul 13 at 18:23