I think there is a bug in ExpectationMatrix.expand.

[JRF-2018]

With ExpectationMatrix.expand now it looks like this (in version 1.10.1):

>>> from sympy.stats import Expectation
>>> from sympy.stats.rv import RandomSymbol, RandomMatrixSymbol

>>> epsilon = RandomSymbol("epsilon")
>>> zeta = RandomMatrixSymbol("zeta", 2, 1)
>>> A = MatrixSymbol("A", 2, 2)

>>> Expectation(A * zeta).expand()
A*ExpectationMatrix(zeta)
>>> Expectation(zeta.T * A).expand()
ExpectationMatrix(zeta.T*A)

To fix this, I think you should modify ExpectationMatrix.expand as follows:

def _ExpectationMatrix_fixed_expand (self, **hints):
  from sympy.stats import Expectation
  from sympy.stats.rv import is_random
  from sympy.core.function import expand as _expand
  expr = self.args[0]
  condition = self._condition
  if not is_random(expr):
      return expr

  if isinstance(expr, Add):
      return Add.fromiter(Expectation(a, condition=condition).expand()
              for a in expr.args)

  expand_expr = _expand(expr)
  if isinstance(expand_expr, Add):
      return Add.fromiter(Expectation(a, condition=condition).expand()
              for a in expand_expr.args)

  elif isinstance(expr, (Mul, MatMul)):
      rv = []
      nonrv = []
      postnon = []

      for a in expr.args:
          if is_random(a):
              if rv:
                  rv.extend(postnon)
              else:
                  nonrv.extend(postnon)
              postnon = []
              rv.append(a)
          elif a.is_Matrix:
              postnon.append(a)
          else:
              nonrv.append(a)

      # In order to avoid infinite-looping (MatMul may call .doit() again),
      # do not rebuild
      if len(nonrv) == 0 and len(postnon) == 0: # <- fixed
          return self
      return Mul.fromiter(nonrv)*Expectation(Mul.fromiter(rv),
              condition=condition)*Mul.fromiter(postnon)

  return self

def fix_ExpectationMatrix_expand ():
  from sympy.stats import ExpectationMatrix
  ExpectationMatrix.expand = _ExpectationMatrix_fixed_expand

>>> fix_ExpectationMatrix_expand()

>>> Expectation(A * zeta).expand()
A*ExpectationMatrix(zeta)
>>> Expectation(zeta.T * A).expand()
ExpectationMatrix(zeta.T)*A

I've only tried it by rewriting a method like above, and I haven't done the whole test, so I chose writing this issue instead of pull request. I also don't want to use pull requests because I'm not familiar with GitHub and I don't speak English well.

[JRF-2018]

I hope ExpectationMatrix can support Transpose as well.

I think the above function should be updated as follows.

def _ExpectationMatrix_fixed_expand (self, **hints):
  from sympy.stats import Expectation
  from sympy.stats.rv import is_random
  from sympy.core.function import expand as _expand
  expr = self.args[0]
  condition = self._condition
  if not is_random(expr):
      return expr

  if isinstance(expr, Add):
      return Add.fromiter(Expectation(a, condition=condition).expand()
              for a in expr.args)

  elif isinstance(expr, Transpose):
      return expr.func(Expectation(expr.args[0],
                                   condition=condition).expand())

  expand_expr = _expand(expr)
  if isinstance(expand_expr, Add):
      return Add.fromiter(Expectation(a, condition=condition).expand()
              for a in expand_expr.args)

  elif isinstance(expand_expr, Transpose):
      return expand_expr.func(Expectation(expand_expr.args[0],
                                          condition=condition).expand())

  elif isinstance(expr, (Mul, MatMul)):
      rv = []
      nonrv = []
      postnon = []

      for a in expr.args:
          if is_random(a):
              if rv:
                  rv.extend(postnon)
              else:
                  nonrv.extend(postnon)
              postnon = []
              rv.append(a)
          elif a.is_Matrix:
              postnon.append(a)
          else:
              nonrv.append(a)

      # In order to avoid infinite-looping (MatMul may call .doit() again),
      # do not rebuild
      if len(nonrv) == 0 and len(postnon) == 0:
          return self
      return Mul.fromiter(nonrv)*Expectation(Mul.fromiter(rv),
              condition=condition).expand()*Mul.fromiter(postnon)

  return self

def fix_ExpectationMatrix_expand ():
  from sympy.stats import ExpectationMatrix
  ExpectationMatrix.expand = _ExpectationMatrix_fixed_expand

>>> fix_ExpectationMatrix_expand()

>>> Expectation(A * zeta).expand()
A*ExpectationMatrix(zeta)
>>> Expectation(zeta.T * A).expand()
ExpectationMatrix(zeta).T*A

[JRF-2018]

I'm sorry, I made a mistake.

The first comment is correct, but the second one was wrong.

I edited the part of the second one from

      if len(nonrv) == 0 and \
         (len(rv) == 0 or len(postnon) == 0):
          return self

to

      if len(nonrv) == 0 and len(postnon) == 0:
          return self