def error_perturbation_single(C, S, R=None):
    r"""Error-perturbation arising from a given sensitivity

    Parameters
    ----------
    C : (M, M) ndarray
        Count matrix
    S : (M, M) ndarray
        Sensitivity matrix
    R : (M, M) ndarray (optional)
        Sensitivity matrix

    Returns
    -------
    var : float
         Variance (covariance) of observable(s)

    """
    cov = tmatrix_cov(C)  # (M, M, M)
    if R is None:
        R = S
    X = S[:, :, np.newaxis] * cov * R[:, np.newaxis, :]
    return X.sum()