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