# Compute the mean, variance and standard deviation of all pairs, triplets and sets of four

My code sample for finding the mean for the case of pairs is:

import numpy as np
Data =[]
with open('testdata.txt') as my_file:
for line in my_file:
Data.append(line)
psize = Data.size*Data.size
p_avg =zeros(psize,float)
for i in range(0,Data.size):
for j in range(0,Data.size):
if i != j:
p_avg[(Data.size*i)+j] = (Data[i]+Data[j])/2.0
else:
p_avg[(Data.size*i)+j] = 0.0
print p_avg


Any comments on this and hints for finding the mean and variance for triplets and sets of four efficiently would be most appreciated.

• Welcome to Code Review, your post looks good, hope you get some good answers! Note that questions about code to be written are in general off topic, but the generalisation from pairs is probably fine, IMO. Oct 5, 2016 at 22:22

If you want to keep all results, this can be done more efficiently using meshgrid. You can create N-dimensional grids of values, flatten them, then take the mean of the result.

So for the 2D case, meshgrid is pretty straightforward:

>>> a = np.arange(5)
>>> [a1, a2] = np.meshgrid(a, a)
>>> print(a1)
[[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]]
>>> print(a2)
[[0 0 0 0 0]
[1 1 1 1 1]
[2 2 2 2 2]
[3 3 3 3 3]
[4 4 4 4 4]]
>>> print(a1.ravel())
[0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4]
>>> print(a2.ravel())
[0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4]


As you can see, the flattened versions of a1 and a2 combined provide every possible combination. This can be extended to higher dimensions just by putting more a values:

>>> a = np.arange(3)
>>> [a1, a2, a3] = np.meshgrid(a, a, a)
>>> b = np.vstack([a1.ravel(), a2.ravel(), a3.ravel()])
>>> print(b)
[[0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2]
[0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2]
[0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2]]


This can be extended to any number of dimensions, allowing of course for memory constraints.

You can generalize the meshgrid function call by using unpacking. These two function calls are equivalent:

np.meshgrid(a, b, c)
np.meshgrid(*[a, b, c])


The * in this case means "unpack the list as function arguments" (also works with tuples, iteratirs, and many other things). You can combine this with the [x]*n syntax, where n is a number. This means "duplicate the contents of list [x] n times". So *4 is similar to [5, 5, 5, 5] ("similar", not "identical", but the differences are not relevant to what we are doing).

So you can combine this to make calling meshgrid with multiple values easier. These three lines of code are the same:

np.meshgrid(a, a, a)
np.meshgrid([a, a, a])
np.meshgrid([a]*3)


You can exclude combinations where all the values are equal using logical indexing. First, you determine whether each element is equal to the value of the first element in that combination. Then you find combinations where that is True for all values. Then you get the columns for which that isn't the case. So something like:

>>> c = b[:, ~(b == b).all(0)]
>>> print(c)
[[0 0 1 1 1 2 2 2 0 0 0 1 1 2 2 2 0 0 0 1 1 1 2 2]
[0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2]
[1 2 0 1 2 0 1 2 0 1 2 0 2 0 1 2 0 1 2 0 1 2 0 1]]


If you want to avoid memory issues you can make use of some internal knowledge of how numpy handles arrays. Arrays are simply continguous (or strided, i.e. continous in steps) blocks of memory. By working with this sort of memory, numpy can make use of very fast libraries for doing calculations on a block of memory at a time. But since this is simply a sequence of bytes, there are different ways to look at this block of memory. The same block could be seen as a 1D array of floats, a 2D array of integers, etc. Numpy can make use of this to avoid making unnecessary copies. So arr.ravel(), for example, "flattens" an array by treating the same underlying data as a 1D array. The catch (or benefit, depending on your use-case) is that changes to a raveled array will affect the original.

meshgrid can make use of this to avoid making unnecessary copies. Consider this:

>>> a = np.arange(10)
>>> b = np.meshgrid(*[a, a, a, a])
>>> c = np.meshgrid(*[a, a, a, a], copy=False)
>>> b.flags['OWNDATA']
True
>>> c.flags['OWNDATA']
False


flags['OWNDATA'] tells you whether the array is a view into the data stored in another array or not. In the case of b, it is, which means each b array takes up memory equal to the size of the array. In the case of c, however, it isn't, which means c takes almost no additional memory, it is still using the same block of memory as a. The catch is that any change to any element to any of the arrays in c will change the value in a and any of the values corresponding to that value in any of the c arrays.

If you combine this with ravel, you can get a view as well, but reshape(-1) is more reliable. Consider this:

>>> c.ravel().flags['OWNDATA']
True
>>> c.reshape(-1).flags['OWNDATA']
False


Slice of arrays, such as a[5:1000] also return views. "Fancy" indexing, where you get multiple elements of an array by index or boolean value, return new arrays. So you can slice out as many values as you want.

Getting the histogram of rows or columns of an array is a bit more difficult. The np.histogram can't do it. What you can do is follow the instructions here to get the unique rows (you need unique columns so just transpose), but when they use the unique function add the return_counts=True argument to also get the counts of each unique row. You will also need to sort each column using d.sort(axis=0) at the beginning so that columns with different orders of the same value are not counted separately.

Note that this doesn't use the highly-optimized numpy array processing libraries, so it will be very slow for large data set.

You can also read all the lines into a numpy array at once using loadtxt

So to implement what you want you can do:

from __future__ import print_function

import numpy as np

n = 4  # This is the number of repetitions
maxn = 10000

vals = np.meshgrid(*[data]*n, copy=False)
vals = np.vstack([val.reshape(-1)[:maxn] for val in vals])
vals = vals[:, ~(vals == vals).all(0)]
p_avg = vals.mean(0)
print(p_avg)


Or as a function:

from __future__ import print_function

import numpy as np

def allmeans(data, n):
vals = np.meshgrid(*[data]*n)
vals = np.vstack([val.ravel() for val in vals])
vals = vals[:, ~(vals == vals).all(0)]
return vals.mean(0)

print(allmeans(data, 4))


Doing variance or standard deviation just requires replacing mean with var and std, respectively. You could even pass it as a function:

from __future__ import print_function

import numpy as np

def allagg(data, n, func):
vals = np.meshgrid(*[data]*n)
vals = np.vstack([val.ravel() for val in vals])
vals = vals[:, ~(vals == vals).all(0)]
return func(vals, 0)

print('Mean: \n', allagg(data, 4, np.mean))
print('Variance: \n', allagg(data, 4, np.var))
print('Standard Deviation: \n', allagg(data, 4, np.std))


NOTE: The differences between [x, x, x] and [x]*3 come into play when x is itself mutable (a data type that can change, like a list or numpy array). Consider this example:

>>> a = [, , ]
>>> b = []*3
>>> a = 5
>>> b = 5
>>> print(a)
[, , ]
>>> print(b)
[, , ]


This difference happens because in a each  is an independent list, while in b each  is the same list. Since you can change lists in-place, changes to a list in one part of your program affect that list everywhere else in the program. So these are more equivalent:

>>> x = 
>>> a = [x, x, x]
>>> b = [x]*3
>>> a = 5
>>> print(a)
[, , ]
>>> print(b)
[, , ]


Notice I didn't even have top change b here, changes to a change every x used everywhere. This can be used to do funny things like create nested infinitely-nested lists:

>>> x = [1, 1]
>>> x = x
>>> x
[[...], 1]


This issue affects lists of numpy arrays as well, but it doesn't matter to use because we aren't changing any of those arrays in-place.

EDIT 1: Change code to exclude duplicated

EDIT 2: Explain unpacking

EDIT 3: Warning about [x]*n syntax

EDIT 4: Added optimization tips and histogram

• Could you explain what np.meshgrid(*[data]*n) is doing?
– Tom
Oct 6, 2016 at 0:06
• @Tom I have edited the answer to include an explanation of that and excluding combinations where all the members are the same. Oct 6, 2016 at 0:16
• Thank you very much for the detailed answer. I have a much better idea of how to use numpy now. Cheers
– Tom
Oct 6, 2016 at 0:17
• @Tom: I added a note at the end about the [x]*n syntax. Please read it, it is a very common source of errors for beginning (and occasionally even experienced) Python programmers. Oct 6, 2016 at 0:25
• Is there a way to truncate the calculation for sets of four to e.g. the first 10,000 such sets? The computer my code has to run on is very old and slow. Also the question asks for histograms of the pairs, triplets and sets of four. I can get a histogram for just the list of numbers easily enough but am not sure how to get the histogram for the generalized cases. Thanks
– Tom
Oct 6, 2016 at 2:27