docs: def solve(a, b):
      """ Find ``x`` such that ``a @ x = b`` for matrix ``a``.
  
      Args:
          a: Two-dimensional, square matrix/array of numbers
              and/or :class:`gvar.GVar`\s.
          b: One-dimensional vector/array of numbers and/or
              :class:`gvar.GVar`\s, or an array of such vectors.
              Requires ``b.shape[0] == a.shape[1]``.
  
      Returns:
          The solution ``x`` of ``a.dot(x) = b``, which is equivalent
          to ``inv(a).dot(b)``.
  
      Raises:
          ValueError: If ``a`` is not square and two-dimensional.
          ValueError: If shape of ``b`` does not match that of ``a``
              (that is ``b.shape[0] != a.shape[1]``).
      """
      amean = gvar.mean(a)
      if amean.ndim != 2 or amean.shape[0] != amean.shape[1]:
          raise ValueError('bad matrix shape: ' + str(a.shape))
      bmean = gvar.mean(b)
      if bmean.shape[0] != a.shape[1]:
          raise ValueError(
              'Mismatch between shapes of a and b: {} {}'.format(a.shape, b.shape)
              )
      # xmean = numpy.linalg.solve(amean, bmean)
      ainv = numpy.linalg.inv(amean)
      xmean = ainv.dot(bmean)
      return xmean + ainv.dot(b-bmean - (a-amean).dot(xmean))