Why using np.mean() and mean() gave me different output number?

Viewed 369

It is interesting to notice that using np.mean() or mean() gave me different outputs.

from statistics import mean
import numpy as np
import matplotlib.pyplot as plt

xs = np.array([1, 2, 3, 4, 5, 6])
ys = np.array([5, 4, 6, 5, 6, 7])

def best_fit_slope(xs, ys):
    numerator = (mean(xs)*mean(ys)) - mean(xs*ys)
    denominator = mean(xs)**2 - mean(xs**2)
    return numerator/denominator 

m = best_fit_slope(xs, ys)
print(m)

output>>> 0.8333333333333334

But if I replace mean() by np.mean() the output >>> 0.42857142857142866.

I followed this video:this video. He just used mean() and gave an output of 0.42857. Can anyone explain why there is a difference? I know most of the linear algebra operations or that involved in arrays, I would prefer to use np.mean().

2 Answers

Interestingly enough, this nuance is not explicitly documented in the formal docs, but it can be deduced from the provided examples.

statistics.mean tries its best to provide output of the same type of the input. When you give it np.array([1, 2, 3, 4, 5, 6]) (array of np.int32) it assumes an int output is expected:

xs = np.array([1, 2, 3, 4, 5, 6])
print(mean(xs))
# 3
print(type(mean(xs)))
# <class 'numpy.int32'>

It is enough to force one value in the array to float to "convince" it we want a float back:

xs = np.array([1.0, 2, 3, 4, 5, 6])
# or np.array([1,2,3,4,5,6],dtype=np.float64) or anyother way that gives `dtype` np.float
print(mean(xs))
# 3.5
print(type(mean(xs)))
# <class 'numpy.float64'>

If we dig deep enough into its implementation we can see where this behaviour comes from. It is using the _sum function which is documented as below:

def _sum(data, start=0):
    """_sum(data [, start]) -> (type, sum, count)

    Return a high-precision sum of the given numeric data as a fraction,
    together with the type to be converted to and the count of items.

    If optional argument ``start`` is given, it is added to the total.
    If ``data`` is empty, ``start`` (defaulting to 0) is returned.


    Examples
    --------

    >>> _sum([3, 2.25, 4.5, -0.5, 1.0], 0.75)
    (<class 'float'>, Fraction(11, 1), 5)

    Some sources of round-off error will be avoided:

    # Built-in sum returns zero.
    >>> _sum([1e50, 1, -1e50] * 1000)
    (<class 'float'>, Fraction(1000, 1), 3000)

    Fractions and Decimals are also supported:

    >>> from fractions import Fraction as F
    >>> _sum([F(2, 3), F(7, 5), F(1, 4), F(5, 6)])
    (<class 'fractions.Fraction'>, Fraction(63, 20), 4)

    >>> from decimal import Decimal as D
    >>> data = [D("0.1375"), D("0.2108"), D("0.3061"), D("0.0419")]
    >>> _sum(data)
    (<class 'decimal.Decimal'>, Fraction(6963, 10000), 4)

    Mixed types are currently treated as an error, except that int is
    allowed.
    """

This is due to how the statistics package tries to give you consistent outputs depending on the numeric type you pass in, so it handles int, float, decimal.Decimal, fractions.Fraction as you would hope. Unfortunately, numpy types don't play well with the python numeric type hierarchy. So we can look at the source code (this is the Python version, your runtime is probably using a fast, C version but they should work equivalently...):

def mean(data):
    """Return the sample arithmetic mean of data.
    >>> mean([1, 2, 3, 4, 4])
    2.8
    >>> from fractions import Fraction as F
    >>> mean([F(3, 7), F(1, 21), F(5, 3), F(1, 3)])
    Fraction(13, 21)
    >>> from decimal import Decimal as D
    >>> mean([D("0.5"), D("0.75"), D("0.625"), D("0.375")])
    Decimal('0.5625')
    If ``data`` is empty, StatisticsError will be raised.
    """
    if iter(data) is data:
        data = list(data)
    n = len(data)
    if n < 1:
        raise StatisticsError('mean requires at least one data point')
    T, total, count = _sum(data)
    assert count == n
    return _convert(total/n, T)

So, essentially it uses a type-aware sum, which returns the type, total, and count. Essentially, the total/count is coerced to T. Notice:

In [28]: T, total, count = statistics._sum(np.array([1,2,3]))

In [29]: T, total, count
Out[29]: (numpy.int64, Fraction(6, 1), 3)

In [30]: total / count
Out[30]: Fraction(2, 1)

In [31]: T(total / count)
Out[31]: 2

Notice, all the objects integers you see here are actually numpy.int64, not vanilla int objects. But why doesn't this happen when we do statistics.mean([1,2,3,4])? Well, because the library was built assuming normal python numeric types, peeking at the _convert function:

def _convert(value, T):
    """Convert value to given numeric type T."""
    if type(value) is T:
        # This covers the cases where T is Fraction, or where value is
        # a NAN or INF (Decimal or float).
        return value
    if issubclass(T, int) and value.denominator != 1:
        T = float
    try:
        # FIXME: what do we do if this overflows?
        return T(value)
    except TypeError:
        if issubclass(T, Decimal):
            return T(value.numerator)/T(value.denominator)
        else:
            raise

You'll notice, it special cases: if issubclass(T, int) and value.denominator != 1, i.e., you have an int and the denominator isn't one, so you want a float:

        T = float

HOWEVER:

In [36]: issubclass(np.int64, int)
Out[36]: False

So, T is just np.int64, and:

In [37]: total / count
Out[37]: Fraction(2, 1)

In [38]: np.int64(total / count)
Out[38]: 2
Related