You can create a 3D array containing your 2D arrays to be averaged, then average along axis=0 using np.mean or np.average (the latter allows for weighted averages):

np.mean( np.array([ old_set, new_set ]), axis=0 )

This averaging scheme can be applied to any (n)-dimensional array, because the created (n+1)-dimensional array will always contain the original arrays to be averaged along its axis=0.

I would use np.mean([ar1,ar2,ar3,ar4...arn], axis=0).

Use the functional form of np.mean:

>>> import numpy as np
>>> arrays = [np.random.random((4,2)) for _ in range(3)]
>>> np.mean(arrays, axis=0)

This converts your list of arrays to a 3D array of shape (3, 4, 2) and then takes the mean along axis 0.

You can also use Python's sum:

>>> sum(arrays)/len(arrays)

For small lists like your example this is actually faster.

Some timings (m is the length of the list):

m: 3   n:4   k: 2
numpy                 0.01291340 ms
python                0.00295936 ms
m: 10   n:100   k: 100
numpy                 0.14189354 ms
python                0.09465128 ms
m: 1000   n:10   k: 10
numpy                 0.43023768 ms
python                0.45201713 ms

Benchmarking code:

import numpy as np

from timeit import timeit
import types

def setup(m, n, k):
    return list(np.random.random((m, n, k)))

def f_numpy(a):
    return np.mean(a, axis=0)

def f_python(a):
    return sum(a)/len(a)

for args in [(3, 4, 2), (10, 100, 100), (1000, 10, 10)]:
    data = setup(*args)
    print('m: {}   n:{}   k: {}'.format(*args))
    for name, func in list(globals().items()):
        if not name.startswith('f_') or not isinstance(func, types.FunctionType):
            continue
        print("{:16s}{:16.8f} ms".format(name[2:], timeit(
            'f(data)', globals={'f':func, 'data':data}, number=1000)))

Your record array from the example above is three dimensional, with shape:

>>> record.shape
(2, 10, 2)

The first dimension corresponds to the 2 iterations of your experiment. To average them, you need to tell np.average to do its thing along axis=0

>>> np.average(record, axis=0)
array([[ 0.        ,  0.45688836],
       [ 0.91377672,  1.37066507],
       [ 1.82755343,  2.28444179],
       [ 2.74133015,  3.19821851],
       [ 3.65510686,  4.11199522],
       [ 4.56888358,  5.02577194],
       [ 5.4826603 ,  5.93954865],
       [ 6.39643701,  6.85332537],
       [ 7.31021373,  7.76710209],
       [ 8.22399044,  8.6808788 ]])

If you know beforehand how many simulations you are going to run, you are better off skipping the list thing altogether and doing something like this:

simulations, sim_rows, sim_cols = 1000000, 10, 2
record = np.empty((simulations, sim_rows, sim_cols))
for j in xrange(simulations) :
    record[j] = np.random.rand(sim_rows, sim_cols)

>>> np.average(record, axis=0)
[[ 0.50021935  0.5000554 ]
 [ 0.50019659  0.50009123]
 [ 0.50008591  0.49973058]
 [ 0.49995812  0.49973941]
 [ 0.49998854  0.49989957]
 [ 0.5002542   0.50027464]
 [ 0.49993122  0.49989623]
 [ 0.50024623  0.49981818]
 [ 0.50005848  0.50016798]
 [ 0.49984452  0.49999112]]

Well, it's easy to do what you describe just by reshaping the array:

import numpy as np
tempList = np.array([[214,295], [215, 294],[229, 226], [229, 227]])
tempList = tempList.reshape( (-1,2,2) )
print(tempList)
print("---")
print( tempList.mean( axis=1 ) )

Output:

[[[214 295]
  [215 294]]

 [[229 226]
  [229 227]]]
---
[[214.5 294.5]
 [229.  226.5]]

To compute some operation slice by slice in numpy, it is very often useful to reshape your array and use extra axes.

To explain the process we'll use here: you can reshape your array, take the mean, reshape it again and take the mean again. Here I assume blocksize is 2

t = np.array([[0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5],[0, 1, 2, 3, 4, 5],[0, 1, 2, 3, 4, 5],[0, 1, 2, 3, 4, 5],[0, 1, 2, 3, 4, 5],])
t = t.reshape([6, 3, 2])
t = np.mean(t, axis=2)
t = t.reshape([3, 2, 3])
np.mean(t, axis=1)

outputs

array([[ 0.5,  2.5,  4.5],
       [ 0.5,  2.5,  4.5],
       [ 0.5,  2.5,  4.5]])

Now that it's clear how this works, you can do it in one pass only:

t = t.reshape([3, 2, 3, 2])
np.mean(t, axis=(1, 3))

works too (and should be quicker since means are computed only once - I guess). I'll let you substitute height/blocksize, width/blocksize and blocksize accordingly.

See @askewcan nice remark on how to generalize this to any dimension.

I often needed this for plotting mean of performance curves with different lengths.

Solved it with simple function (based on answer of @unutbu):

def tolerant_mean(arrs):
    lens = [len(i) for i in arrs]
    arr = np.ma.empty((np.max(lens),len(arrs)))
    arr.mask = True
    for idx, l in enumerate(arrs):
        arr[:len(l),idx] = l
    return arr.mean(axis = -1), arr.std(axis=-1)

y, error = tolerant_mean(list_of_ys_diff_len)
ax.plot(np.arange(len(y))+1, y, color='green')

So applying that function to the list of above-plotted curves yields the following: