3
\$\begingroup\$

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]))
\$\endgroup\$
1
\$\begingroup\$

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}$$

\$\endgroup\$
  • 1
    \$\begingroup\$ 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); \$\endgroup\$ – Erel Segal-Halevi Jul 13 at 19:01
  • \$\begingroup\$ 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. \$\endgroup\$ – dfhwze Jul 13 at 19:07
  • \$\begingroup\$ @ErelSegal-Halevi Btw, read the answers on this question for a single pass calculation of both mean and variance: codereview.stackexchange.com/questions/185450/… \$\endgroup\$ – dfhwze Jul 13 at 19:41
1
\$\begingroup\$

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)
\$\endgroup\$
  • 1
    \$\begingroup\$ I think whatever_that_is_case is called snake_case (en.wikipedia.org/wiki/Snake_case). \$\endgroup\$ – AlexV Jul 12 at 16:33
  • \$\begingroup\$ @AlexV thanks that's fixed \$\endgroup\$ – IEatBagels Jul 12 at 16:34
  • 1
    \$\begingroup\$ It seems like the Python style guide itself refers to that style as lower_case_with_underscores. But it doesn't really matter :-D \$\endgroup\$ – AlexV Jul 12 at 16:46
  • \$\begingroup\$ 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? \$\endgroup\$ – Erel Segal-Halevi Jul 13 at 18:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.