docs: def get_covariance(datargs, outargs, vargs, datvar, outvar):
          """
          Get covariance matrix.
  
          :param datargs: data arguments
          :param outargs: output arguments
          :param vargs: variable arguments
          :param datvar: variance of data arguments
          :param outvar: variance of output arguments
          :return: covariance
          """
          # number of formula arguments that are not constant
          argn = len(vargs)
          # number of observations must be the same for all vargs
          nobs = 1
          for m in xrange(argn):
              a = vargs[m]
              try:
                  a = datargs[a]
              except (KeyError, TypeError):
                  a = outargs[a]
                  avar = outvar[a]
              else:
                  avar = datvar[a]
              for n in xrange(argn):
                  b = vargs[n]
                  try:
                      b = datargs[b]
                  except (KeyError, TypeError):
                      b = outargs[b]
                  c = avar.get(b, 0.0)
                  try:
                      nobs = max(nobs, len(c))
                  except (TypeError, ValueError):
                      LOGGER.debug('c of %s vs %s = %g', a, b, c)
          # covariance matrix is initially zeros
          cov = np.zeros((nobs, argn, argn))
          # loop over arguments in both directions, fill in covariance
          for m in xrange(argn):
              a = vargs[m]
              try:
                  a = datargs[a]
              except (KeyError, TypeError):
                  a = outargs[a]
                  avar = outvar[a]
              else:
                  avar = datvar[a]
              for n in xrange(argn):
                  b = vargs[n]
                  try:
                      b = datargs[b]
                  except (KeyError, TypeError):
                      b = outargs[b]
                  cov[:, m, n] = avar.get(b, 0.0)
          if nobs == 1:
              cov = cov.squeeze()  # squeeze out any extra dimensions
          LOGGER.debug('covariance:\n%r', cov)
          return cov