Poor numerical consistency with numpy for reductions with dtype=`float16` specified

Reductions like COO.sum and COO.mean fail to match numpy in all cases when dtype=numpy.float16 is specified. For example:

import sparse
import numpy as np
x = np.array([[[    0,     0,  4526,     0],
               [    0,     0,   -37,     0],
               [ 8372,     0,  7915,     0]],
              [[    0,     0,     0,     0],
               [    0,     0, -7917,     0],
               [-9719,     0,     0,     0]]], dtype='i4')

s = sparse.COO.from_numpy(x)
res = s.sum(axis=(0, 2), dtype='f2')
sol = x.sum(axis=(0, 2), dtype='f2')
print(res.todense())
print(sol)

outputs:

[ 4530. -7950.  6564.]
[ 4530. -7950.  6570.]

It's not clear if each of these will need to be fixed per-method, or if there's a general fix in the reduce code. For sum, numpy interprets x.sum(dtype='f2') the same as x.astype('f2').sum(), but this is not true for x.mean(dtype='f2').

Start at the COO.sum and COO.mean reduction entry points and reproduce the provided float16 examples against NumPy. Trace how the dtype argument is handled for reductions over multiple axes, then verify that sparse results match the corresponding NumPy results for the shown cases.