code: def combine(self, other, mean_free=False):
          
          w1 = self.w
          w2 = other.w
          w = w1 + w2
          # TODO: fix this div by zero error
          q = w2 / w1
          dsx = q * self.sx - other.sx
          dsy = q * self.sy - other.sy
          # update
          self.w = w1 + w2
          self.sx = self.sx + other.sx
          self.sy = self.sy + other.sy
          #
          if mean_free:
              if len(self.Mxy.shape) == 1:  # diagonal only
                  d = dsx*dsy
              else:
                  d = np.outer(dsx, dsy)
              self.Mxy += other.Mxy + (w1 / (w2 * w)) * d
          else:
              self.Mxy += other.Mxy
          return self
docs: References
          ----------
          [1] http://i.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf