Repository: statsmodels/statsmodels (Python)

Issue:
BUG: RLM fit with start_params raises if only one parameter
if `start_params = np.array([5])` one element, 1-dim array, then the squeeze removes the 1-dim which then raises exception in shape check

```
statsmodels_gh\statsmodels\statsmodels\robust\robust_linear_model.py in fit(self, maxiter, tol, scale_est, init, cov, update_scale, conv, start_params, start_scale)
    264         else:
    265             start_params = np.asarray(start_params, dtype=np.double).squeeze()
--> 266             if (start_params.shape[0] != self.exog.shape[1] or
    267                     start_params.ndim != 1):
    268                 raise ValueError('start_params must by a 1-d array with {} '

IndexError: tuple index out of range
```

I ran into this for the constant only regression for resistant estimators in #9227

I don't know why there is the squeeze.
If we want or need to keep the squeeze, then a fix is to add
`start_params = np.atleast_1d(start_params),`

I will fix it this way in my PR, where I need it now


Hints:

Patch:
diff --git a/statsmodels/robust/covariance.py b/statsmodels/robust/covariance.py
--- a/statsmodels/robust/covariance.py
+++ b/statsmodels/robust/covariance.py
@@ -28,9 +28,15 @@
 from scipy import stats, linalg
 from scipy.linalg.lapack import dtrtri
 from .scale import mad
+import statsmodels.robust.norms as rnorms
+import statsmodels.robust.scale as rscale
+import statsmodels.robust.tools as rtools
 from statsmodels.tools.testing import Holder
 
+from statsmodels.stats.covariance import corr_rank, corr_normal_scores
+
 mad0 = lambda x: mad(x, center=0)  # noqa: E731
+median = lambda x: np.median(x, axis=0)  # noqa: E731
 
 
 # from scikit-learn
@@ -311,6 +317,7 @@ def mahalanobis(data, cov=None, cov_inv=None, sqrt=False):
     """Mahalanobis distance squared
 
     Note: this is without taking the square root.
+    assumes data is already centered.
 
     Parameters
     ----------
@@ -405,13 +412,15 @@ def cov_gk(data, scale_func=mad):
 def cov_ogk(data, maxiter=2, scale_func=mad, cov_func=cov_gk,
             loc_func=lambda x: np.median(x, axis=0), reweight=0.9,
             rescale=True, rescale_raw=True, ddof=1):
-    """orthogonalized Gnanadesikan and Kettenring covariance estimator
+    """Orthogonalized Gnanadesikan and Kettenring covariance estimator.
 
     Based on Maronna and Zamar 2002
 
     Parameters
     ----------
-    data : array_like, 2-D
+    data : array-like
+        Multivariate data set with observation in rows and variables in
+        columns.
     maxiter : int
         Number of iteration steps. According to Maronna and Zamar the
         estimate doesn't improve much after the second iteration and the
@@ -590,12 +599,10 @@ def cov_tyler(data, start_cov=None, normalize=False, maxiter=100, eps=1e-13):
     ----------
     .. [1] Tyler, David E. “A Distribution-Free M-Estimator of Multivariate
        Scatter.” The Annals of Statistics 15, no. 1 (March 1, 1987): 234–51.
-
     .. [2] Soloveychik, I., and A. Wiesel. 2014. Tyler's Covariance Matrix
        Estimator in Elliptical Models With Convex Structure.
        IEEE Transactions on Signal Processing 62 (20): 5251-59.
        doi:10.1109/TSP.2014.2348951.
-
     .. [3] Ollila, Esa, Daniel P. Palomar, and Frederic Pascal.
        “Affine Equivariant Tyler’s M-Estimator Applied to Tail Parameter
        Learning of Elliptical Distributions.” arXiv, May 7, 2023.
@@ -1083,7 +1090,7 @@ def _cov_iter(data, weights_func, weights_args=None, cov_init=None,
     return res
 
 
-def _cov_starting(data, standardize=False, quantile=0.5):
+def _cov_starting(data, standardize=False, quantile=0.5, retransform=False):
     """compute some robust starting covariances
 
     The returned covariance matrices are intended as starting values
@@ -1127,7 +1134,7 @@ def _cov_starting(data, standardize=False, quantile=0.5):
 
     cov_all = []
     d = mahalanobis(xs, cov=None, cov_inv=np.eye(k_vars))
-    percentiles = [(k_vars+2) / nobs * 100 * 2, 25, 50]
+    percentiles = [(k_vars+2) / nobs * 100 * 2, 25, 50, 85]
     cutoffs = np.percentile(d, percentiles)
     for p, cutoff in zip(percentiles, cutoffs):
         xsp = xs[d < cutoff]
@@ -1149,15 +1156,929 @@ def _cov_starting(data, standardize=False, quantile=0.5):
         c03 = _cov_iter(xs, weights_quantile, weights_args=(quantile,),
                         rescale="med", cov_init=c02.cov, maxiter=100)
 
-        if standardize:
+        if not standardize or not retransform:
             cov_all.extend([c0, c01, c02, c03])
         else:
             # compensate for initial rescaling
+            # TODO: this does not return list of Holder anymore
             s = np.outer(std, std)
             cov_all.extend([r.cov * s for r in [c0, c01, c02, c03]])
 
     c2 = cov_ogk(xs)
     cov_all.append(c2)
 
+    c2raw = Holder(
+            cov=c2.cov_raw,
+            mean=c2.loc_raw * std + center,
+            method="ogk_raw",
+            )
+    cov_all.append(c2raw)
+
+    z_tanh = np.tanh(xs)
+    c_th = Holder(
+            cov=np.corrcoef(z_tanh.T),  # not consistently scaled for cov
+            mean=center,  # TODO: do we add inverted mean z_tanh ?
+            method="tanh",
+            )
+    cov_all.append(c_th)
+
+    x_spatial = xs / np.sqrt(np.sum(xs**2, axis=1))[:, None]
+    c_th = Holder(
+            cov=np.cov(x_spatial.T),
+            mean=center,
+            method="spatial",
+            )
+    cov_all.append(c_th)
+
+    c_th = Holder(
+            # not consistently scaled for cov
+            # cov=stats.spearmanr(xs)[0], # not correct shape if k=1 or 2
+            cov=corr_rank(xs),  # always returns matrix, np.corrcoef result
+            mean=center,
+            method="spearman",
+            )
+    cov_all.append(c_th)
+
+    c_ns = Holder(
+            cov=corr_normal_scores(xs),  # not consistently scaled for cov
+            mean=center,  # TODO: do we add inverted mean z_tanh ?
+            method="normal-scores",
+            )
+    cov_all.append(c_ns)
+
     # TODO: rescale back to original space using center and std
     return cov_all
+
+
+# ####### Det, CovDet and helper functions, might be moved to separate module
+
+
+class _Standardize():
+    """Robust standardization of random variable
+
+    """
+
+    def __init__(self, x, func_center=None, func_scale=None):
+        # naming mean or center
+        # maybe also allow str func_scale for robust.scale, e.g. for not
+        #    vectorized Qn
+        if func_center is None:
+            center = np.median(x, axis=0)
+        else:
+            center = func_center(x)
+        xdm = x - center
+        if func_scale is None:
+            scale = mad(xdm, center=0)
+        else:
+            # assumes vectorized
+            scale = func_scale(x)
+
+        self.x_stand = xdm / scale
+
+        self.center = center
+        self.scale = scale
+
+    def transform(self, x):
+        return (x - self.center) / self.scale
+
+    def untransform_mom(self, m, c):
+        mean = self.center + m
+        cov = c * np.outer(self.scale, self.scale)
+        return mean, cov
+
+
+def _orthogonalize_det(x, corr, loc_func, scale_func):
+    """Orthogonalize
+
+    This is a simplified version of the OGK method.
+    version from DetMCD works on zscored data
+    (does not return mean and cov of original data)
+    so we drop the compensation for scaling in zscoring
+
+    z is the data here, zscored with robust estimators,
+    e.g. median and Qn in DetMCD
+    """
+    evals, evecs = np.linalg.eigh(corr)  # noqa: F841
+    z = x.dot(evecs)
+    transf0 = evecs
+
+    scale_z = scale_func(z)  # scale of principal components
+    cov = (transf0 * scale_z**2).dot(transf0.T)
+    # extra step in DetMCD, sphering data with new cov to compute center
+    # I think this is equivalent to scaling z
+    # loc_z = loc_func(z / scale_z) * scale_z  # center of principal components
+    # loc = (transf0 * scale_z).dot(loc_z)
+    transf1 = (transf0 * scale_z).dot(transf0.T)
+    # transf1inv = (transf0 * scale_z**(-1)).dot(transf0.T)
+
+    # loc = loc_func(x @ transf1inv) @ transf1
+    loc = loc_func((z / scale_z).dot(transf0.T)) @ transf1
+
+    return loc, cov
+
+
+def _get_detcov_startidx(z, h, options_start=None, methods_cov="all"):
+    """Starting sets for deterministic robust covariance estimators.
+
+    These are intended as starting sets for DetMCD, DetS and DetMM.
+    """
+
+    if options_start is None:
+        options_start = {}
+
+    loc_func = options_start.get("loc_func", median)
+    scale_func = options_start.get("scale_func", mad)
+    z = (z - loc_func(z)) / scale_func(z)
+
+    if np.squeeze(z).ndim == 1:
+        # only one random variable
+        z = np.squeeze(z)
+        nobs = z.shape[0]
+        idx_sel = np.argpartition(np.abs(z), h)[:h]
+        idx_all = [(idx_sel, "abs-resid")]
+        # next uses symmetric equal-tail trimming
+        idx_sorted = np.argsort(z)
+        h_tail = (nobs - h) // 2
+        idx_all.append((idx_sorted[h_tail : h_tail + h], "trimmed-tail"))
+        return idx_all
+
+    # continue if more than 1 random variable
+    cov_all = _cov_starting(z, standardize=False, quantile=0.5)
+
+    # orthogonalization step
+    idx_all = []
+    for c in cov_all:
+        if not hasattr(c, "method"):
+            continue
+        method = c.method
+        mean, cov = _orthogonalize_det(z, c.cov, loc_func, scale_func)
+        d = mahalanobis(z, mean, cov)
+        idx_sel = np.argpartition(d, h)[:h]
+        idx_all.append((idx_sel, method))
+
+    return idx_all
+
+
+class CovM:
+    """M-estimator for multivariate Mean and Scatter.
+
+    Interface incomplete and experimental.
+
+    Parameters
+    ----------
+    data : array-like
+        Multivariate data set with observation in rows and variables in
+        columns.
+    norm_mean : norm instance
+        If None, then TukeyBiweight norm is used.
+        (Currently no other norms are supported for calling the initial
+        S-estimator)
+    norm_scatter : None or norm instance
+        If norm_scatter is None, then the norm_mean will be used.
+    breakdown_point : float in (0, 0.5]
+        Breakdown point for first stage S-estimator.
+    scale_bias : None or float
+        Must currently be provided if norm_mean is not None.
+    method : str
+        Currently only S-estimator has automatic selection of scale function.
+    """
+
+    def __init__(self, data, norm_mean=None, norm_scatter=None,
+                 scale_bias=None, method="S"):
+        # todo: method defines how norm_mean and norm_scatter are linked
+        #       currently I try for S-estimator
+
+        if method.lower() not in ["s"]:
+            msg = f"method {method} option not recognize or implemented"
+            raise ValueError(msg)
+
+        self.data = np.asarray(data)
+        self.k_vars = k_vars = self.data.shape[1]
+
+        # Todo: check interface for scale bias
+        self.scale_bias = scale_bias
+        if norm_mean is None:
+            norm_mean = rnorms.TukeyBiweight()
+            c = rtools.tuning_s_cov(norm_mean, k_vars, breakdown_point=0.5)
+            norm_mean._set_tuning_param(c, inplace=True)
+            self.scale_bias = rtools.scale_bias_cov_biw(c, k_vars)[0]
+
+        self.norm_mean = norm_mean
+        if norm_scatter is None:
+            self.norm_scatter = self.norm_mean
+        else:
+            self.norm_scatter = norm_scatter
+
+        self.weights_mean = self.norm_mean.weights
+        self.weights_scatter = self.weights_mean
+        # self.weights_scatter = lambda d: self.norm_mean.rho(d) / d**2
+        # this is for S-estimator, M-scale
+        self.rho = self.norm_scatter.rho
+
+    def _fit_mean_shape(self, mean, shape, scale):
+        """Estimate mean and shape in iteration step.
+
+        This does only one step.
+
+        Parameters
+        ----------
+        mean : ndarray
+            Starting value for mean
+        shape : ndarray
+            Starting value for shape matrix.
+        scale : float
+            Starting value for scale.
+
+        Returns
+        -------
+        Holder instance with updated estimates.
+        """
+        d = mahalanobis(self.data - mean, shape, sqrt=True) / scale
+        weights_mean = self.weights_mean(d)
+        weights_cov = self.weights_scatter(d)
+
+        res = cov_weighted(
+            self.data,
+            weights=weights_mean,
+            center=None,
+            weights_cov=weights_cov,
+            weights_cov_denom="det",
+            ddof=1,
+            )
+        return res
+
+    def _fit_scale(self, maha, start_scale=None, maxiter=100, rtol=1e-5,
+                   atol=1e-5):
+        """Estimate iterated M-scale.
+
+        Parameters
+        ----------
+        maha : ndarray
+        start_scale : None or float
+            Starting scale. If it is None, the mad of maha wi
+        maxiter : int
+            Maximum iterations to compute M-scale
+        rtol, atol : float
+            Relative and absolute convergence criteria for scale used with
+            allclose.
+
+        Returns
+        -------
+        float : scale estimate
+        """
+        if start_scale is None:
+            # TODO: this does not really make sense
+            # better scale to median of maha and chi or chi2
+            start_scale = mad(maha)
+
+        scale = rscale._scale_iter(
+            maha,
+            scale0=start_scale,
+            maxiter=maxiter,
+            rtol=rtol,
+            atol=atol,
+            meef_scale=self.rho,
+            scale_bias=self.scale_bias,
+            )
+        return scale
+
+    def fit(self, start_mean=None, start_shape=None, start_scale=None,
+            maxiter=100, update_scale=True):
+        """Estimate mean, shape and scale parameters with MM-estimator.
+
+        Parameters
+        ----------
+        start_mean : None or float
+            Starting value for mean, center.
+            If None, then median is used.
+        start_shape : None or 2-dim ndarray
+            Starting value of shape matrix, i.e. scatter matrix normalized
+            to det(scatter) = 1.
+            If None, then scaled covariance matrix of data is used.
+        start_scale : None or float.
+            Starting value of scale.
+        maxiter : int
+            Maximum number of iterations.
+        update_scale : bool
+            If update_scale is False, then
+
+        Returns
+        -------
+        results instance with mean, shape, scale, cov and other attributes.
+
+        Notes
+        -----
+        If start_scale is provided and update_scale is False, then this is
+        an M-estimator with a predetermined scale as used in the second
+        stage of an MM-estimator.
+
+        """
+
+        converged = False
+
+        if start_scale is not None:
+            scale_old = start_scale
+        else:
+            scale_old = 1
+            # will be reset if start_shape is also None.
+        if start_mean is not None:
+            mean_old = start_mean
+        else:
+            mean_old = np.median(self.data, axis=0)
+        if start_shape is not None:
+            shape_old = start_shape
+        else:
+            shape_old = np.cov(self.data.T)
+            scale = np.linalg.det(shape_old) ** (1 / self.k_vars)
+            shape_old /= scale
+            if start_scale is not None:
+                scale_old = scale
+
+        if update_scale is False:
+            scale = start_scale
+
+        for i in range(maxiter):
+            shape, mean = self._fit_mean_shape(mean_old, shape_old, scale_old)
+            d = mahalanobis(self.data - mean, shape, sqrt=True)
+            if update_scale:
+                scale = self._fit_scale(d, start_scale=scale_old, maxiter=10)
+
+            if (np.allclose(scale, scale_old, rtol=1e-5) and
+                    np.allclose(mean, mean_old, rtol=1e-5) and
+                    np.allclose(shape, shape_old, rtol=1e-5)
+                    ):  # noqa E124
+                converged = True
+                break
+            scale_old = scale
+            mean_old = mean
+            shape_old = shape
+
+        maha = mahalanobis(self.data - mean, shape / scale, sqrt=True)
+
+        res = Holder(
+            mean=mean,
+            shape=shape,
+            scale=scale,
+            cov=shape * scale**2,
+            converged=converged,
+            n_iter=i,
+            mahalanobis=maha,
+            )
+        return res
+
+
+class CovDetMCD:
+    """Minimum covariance determinant estimator with deterministic starts.
+
+    preliminary version
+
+    reproducability:
+
+    This uses deterministic starting sets and there is no randomness in the
+    estimator.
+    However, this will not be reprodusible across statsmodels versions
+    when the methods for starting sets or tuning parameters for the
+    optimization change.
+
+    Parameters
+    ----------
+    data : array-like
+        Multivariate data set with observation in rows and variables in
+        columns.
+
+    Notes
+    -----
+    The correction to the scale to take account of trimming in the reweighting
+    estimator is based on the chisquare tail probability.
+    This differs from CovMcd in R which uses the observed fraction of
+    observations above the metric trimming threshold.
+
+    References
+    ----------
+    ..[1] Hubert, Mia, Peter Rousseeuw, Dina Vanpaemel, and Tim Verdonck. 2015.
+       “The DetS and DetMM Estimators for Multivariate Location and Scatter.”
+       Computational Statistics & Data Analysis 81 (January): 64–75.
+       https://doi.org/10.1016/j.csda.2014.07.013.
+    ..[2] Hubert, Mia, Peter J. Rousseeuw, and Tim Verdonck. 2012. “A
+       Deterministic Algorithm for Robust Location and Scatter.” Journal of
+       Computational and Graphical Statistics 21 (3): 618–37.
+       https://doi.org/10.1080/10618600.2012.672100.
+    """
+
+    def __init__(self, data):
+        # no options yet, methods were written as functions
+        self.data = np.asarray(data)
+
+    def _cstep(self, x, mean, cov, h, maxiter=2, tol=1e-8):
+        """C-step for mcd iteration
+
+        x is data, perc is percentile h / nobs, don't need perc when we
+        use np.argpartition
+        requires starting mean and cov
+        """
+
+        converged = False
+
+        for _ in range(maxiter):
+            d = mahalanobis(x - mean, cov)
+            idx_sel = np.argpartition(d, h)[:h]
+            x_sel = x[idx_sel]
+            mean = x_sel.mean(0)
+            cov_new = np.cov(x_sel.T, ddof=1)
+
+            if ((cov - cov_new)**2).mean() < tol:
+                cov = cov_new
+                converged = True
+                break
+
+            cov = cov_new
+
+        return mean, cov, converged
+
+    def _fit_one(self, x, idx, h, maxiter=2, mean=None, cov=None):
+        """Compute mcd for one starting set of observations.
+
+        Parameters
+        ----------
+        x : ndarray
+            Data.
+        idx : ndarray
+            Indices or mask of observation in starting set, used as ``x[idx]``
+        h : int
+            Number of observations in evaluation set for cov.
+        maxiter : int
+            Maximum number of c-steps.
+
+        Returns
+        -------
+        mean : ndarray
+            Estimated mean.
+        cov : ndarray
+            Estimated covariance.
+        det : float
+            Determinant of estimated covariance matrix.
+
+        Notes
+        -----
+        This does not do any preprocessing of the data and returns the
+        empirical mean and covariance of evaluation set of the data ``x``.
+        """
+        if idx is not None:
+            x_sel = x[idx]
+        else:
+            x_sel = x
+        if mean is None:
+            mean = x_sel.mean(0)
+        if cov is None:
+            cov = np.cov(x_sel.T, ddof=1)
+
+        # updated with c-step
+        mean, cov, conv = self._cstep(x, mean, cov, h, maxiter=maxiter)
+        det = np.linalg.det(cov)
+
+        return mean, cov, det, conv
+
+    def fit(self, h, *, h_start=None, mean_func=None, scale_func=None,
+            maxiter=100, options_start=None, reweight=True,
+            trim_frac=0.975, maxiter_step=100):
+        """
+        Compute minimum covariance determinant estimate of mean and covariance.
+
+        x : array-like
+            Data with observation in rows and variables in columns.
+        h : int
+            Number of observations in evaluation set for minimimizing
+            determinant.
+        h_start : int
+            Number of observations used in starting mean and covariance.
+        mean_func, scale_func : callable or None.
+            Mean and scale function for initial standardization.
+            Current defaults, if they are None, are median and mad, but
+            default scale_func will likely change.
+        options_start : None or dict
+            Options for the starting estimators.
+            currently not used
+            TODO: which options? e.g. for OGK
+        reweight : bool
+            If reweight is true, then a reweighted estimator is returned. The
+            reweighting is based on a chisquare trimming of Mahalanobis
+            distances. The raw results are in the ``results_raw`` attribute.
+        trim_frac : float in (0, 1)
+            Trim fraction used if reweight is true. Used to compute quantile
+            of chisquare distribution with tail probability 1 - trim_frac.
+        maxiter_step : int
+            Number of iteration in the c-step.
+            In the current implementation a small maxiter in the c-step does
+            not find the optimal solution.
+
+        Returns
+        -------
+        Holder instance with results
+        """
+
+        x = self.data
+        nobs, k_vars = x.shape
+
+        if h is None:
+            h = (nobs + k_vars + 1) // 2  # check with literature
+        if mean_func is None:
+            mean_func = lambda x: np.median(x, axis=0)  # noqa
+        if scale_func is None:
+            scale_func = mad
+        if options_start is None:
+            options_start = {}
+        if h_start is None:
+            nobs, k_vars = x.shape
+            h_start = max(nobs // 2 + 1, k_vars + 1)
+
+        m = mean_func(x)
+        s = scale_func(x)
+        z = (x - m) / s
+        # get initial mean, cov of standardized data, we only need ranking
+        # of obs
+        starts = _get_detcov_startidx(z, h_start, options_start)
+
+        fac_trunc = coef_normalize_cov_truncated(h / nobs, k_vars)
+
+        res = {}
+        for ii, ini in enumerate(starts):
+            idx_sel, method = ini
+            mean, cov, det, _ = self._fit_one(x, idx_sel, h,
+                                              maxiter=maxiter_step)
+            res[ii] = Holder(
+                mean=mean,
+                cov=cov * fac_trunc,
+                det_subset=det,
+                method=method,
+                )
+
+        det_all = np.array([i.det_subset for i in res.values()])
+        idx_best = np.argmin(det_all)
+        best = res[idx_best]
+        # mean = best.mean
+        # cov = best.cov
+
+        # need to c-step to convergence for best,
+        # is with best 2 in original DetMCD
+        if maxiter_step < maxiter:
+            mean, cov, det, conv = self._fit_one(x, None, h, maxiter=maxiter,
+                                                 mean=best.mean, cov=best.cov)
+            best = Holder(
+                mean=mean,
+                cov=cov * fac_trunc,
+                det_subset=det,
+                method=method,
+                converged=conv,
+                )
+
+        # include extra info in returned Holder instance
+        best.det_all = det_all
+        best.idx_best = idx_best
+
+        best.tmean = m
+        best.tscale = s
+
+        if reweight:
+            cov, mean = _reweight(x, best.mean, best.cov, trim_frac=trim_frac,
+                                  ddof=1)
+            fac_trunc = coef_normalize_cov_truncated(trim_frac, k_vars)
+            best_w = Holder(
+                mean=mean,
+                cov=cov * fac_trunc,
+                # det_subset=det,
+                method=method,
+                results_raw=best,
+                )
+
+            return best_w
+        else:
+            return best  # is Holder instance already
+
+
+class CovDetS:
+    """S-estimator for mean and covariance with deterministic starts.
+
+    Parameters
+    ----------
+    data : array-like
+        Multivariate data set with observation in rows and variables in
+        columns.
+    norm : norm instance
+        If None, then TukeyBiweight norm is used.
+        (Currently no other norms are supported for calling the initial
+        S-estimator)
+    breakdown_point : float in (0, 0.5]
+        Breakdown point for first stage S-estimator.
+
+    Notes
+    -----
+    Reproducability:
+
+    This uses deterministic starting sets and there is no randomness in the
+    estimator.
+    However, the estimates may not be reproducable across statsmodels versions
+    when the methods for starting sets or default tuning parameters for the
+    optimization change. With different starting sets, the estimate can
+    converge to a different local optimum.
+
+    References
+    ----------
+    ..[1] Hubert, Mia, Peter Rousseeuw, Dina Vanpaemel, and Tim Verdonck. 2015.
+       “The DetS and DetMM Estimators for Multivariate Location and Scatter.”
+       Computational Statistics & Data Analysis 81 (January): 64–75.
+       https://doi.org/10.1016/j.csda.2014.07.013.
+    ..[2] Hubert, Mia, Peter J. Rousseeuw, and Tim Verdonck. 2012. “A
+       Deterministic Algorithm for Robust Location and Scatter.” Journal of
+       Computational and Graphical Statistics 21 (3): 618–37.
+       https://doi.org/10.1080/10618600.2012.672100.
+    """
+
+    def __init__(self, data, norm=None, breakdown_point=0.5):
+        # no options yet, methods were written as functions
+        self.data = np.asarray(data)
+        self.nobs, k_vars = self.data.shape
+        self.k_vars = k_vars
+
+        # self.scale_bias = scale_bias
+        if norm is None:
+            norm = rnorms.TukeyBiweight()
+            c = rtools.tuning_s_cov(norm, k_vars, breakdown_point=0.5)
+            norm._set_tuning_param(c, inplace=True)
+            self.scale_bias = rtools.scale_bias_cov_biw(c, k_vars)[0]
+        else:
+            raise NotImplementedError("only Biweight norm is supported")
+
+        self.norm = norm
+
+        self.mod = CovM(data, norm_mean=norm, norm_scatter=norm,
+                        scale_bias=self.scale_bias, method="S")
+
+    def _get_start_params(self, idx):
+        """Starting parameters from a subsample given by index
+
+        Parameters
+        ----------
+        idx : ndarray
+            Index used to select observations from the data. The index is used
+            for numpy arrays, so it can be either a boolean mask or integers.
+
+        Returns
+        -------
+        mean : ndarray
+            Mean of subsample
+        shape : ndarray
+            The shape matrix of the subsample which is the covariance
+            normalized so that determinant of shape is one.
+        scale : float
+            Scale of subsample, computed so that cov = shape * scale.
+        """
+        x_sel = self.data[idx]
+        k = x_sel.shape[1]
+
+        mean = x_sel.mean(0)
+        cov = np.cov(x_sel.T)
+
+        scale2 = np.linalg.det(cov) ** (1 / k)
+        shape = cov / scale2
+        scale = np.sqrt(scale2)
+        return mean, shape, scale
+
+    def _fit_one(self, mean=None, shape=None, scale=None, maxiter=100):
+        """Compute local M-estimator for one starting set of observations.
+
+        Parameters
+        ----------
+        x : ndarray
+            Data.
+        idx : ndarray
+            Indices or mask of observation in starting set, used as ``x[idx]``
+        h : int
+            Number of observations in evaluation set for cov.
+        maxiter : int
+            Maximum number of c-steps.
+
+        Returns
+        -------
+        mean : ndarray
+            Estimated mean.
+        cov : ndarray
+            Estimated covariance.
+        det : float
+            Determinant of estimated covariance matrix.
+
+        Notes
+        -----
+        This uses CovM to solve for the local optimum for given starting
+        values.
+        """
+
+        res = self.mod.fit(
+            start_mean=mean,
+            start_shape=shape,
+            start_scale=scale,
+            maxiter=maxiter,
+            update_scale=True
+            )
+
+        return res
+
+    def fit(self, *, h_start=None, mean_func=None, scale_func=None,
+            maxiter=100, options_start=None, maxiter_step=5):
+        """Compute S-estimator of mean and covariance.
+
+        Parameters
+        ----------
+        h_start : int
+            Number of observations used in starting mean and covariance.
+        mean_func, scale_func : callable or None.
+            Mean and scale function for initial standardization.
+            Current defaults, if they are None, are median and mad, but
+            default scale_func will likely change.
+        options_start : None or dict
+           Options for the starting estimators.
+           TODO: which options? e.g. for OGK
+
+        Returns
+        -------
+        Holder instance with results
+        """
+
+        x = self.data
+        nobs, k_vars = x.shape
+
+        if mean_func is None:
+            mean_func = lambda x: np.median(x, axis=0)  # noqa
+        if scale_func is None:
+            scale_func = mad
+        if options_start is None:
+            options_start = {}
+        if h_start is None:
+            nobs, k_vars = x.shape
+            h_start = max(nobs // 2 + 1, k_vars + 1)
+
+        m = mean_func(x)
+        s = scale_func(x)
+        z = (x - m) / s
+        # get initial mean, cov of standardized data, we only need ranking
+        # of obs
+        starts = _get_detcov_startidx(z, h_start, options_start)
+
+        res = {}
+        for ii, ini in enumerate(starts):
+            idx_sel, method = ini
+            mean0, shape0, scale0 = self._get_start_params(idx_sel)
+            res_i = self._fit_one(
+                mean=mean0,
+                shape=shape0,
+                scale=scale0,
+                maxiter=maxiter_step,
+                )
+
+            res_i.method = method
+            res[ii] = res_i
+
+        scale_all = np.array([i.scale for i in res.values()])
+        idx_best = np.argmin(scale_all)
+        best = res[idx_best]
+        # mean = best.mean
+        # cov = best.cov
+
+        # need to c-step to convergence for best,
+        # is with best 2 in original DetMCD
+        if maxiter_step < maxiter:
+            best = self._fit_one(
+                mean=best.mean,
+                shape=best.shape,
+                scale=best.scale,
+                maxiter=maxiter,
+                )
+
+        # include extra info in returned Holder instance
+        best.scale_all = scale_all
+        best.idx_best = idx_best
+
+        best.tmean = m
+        best.tscale = s
+
+        return best  # is Holder instance already
+
+
+class CovDetMM():
+    """MM estimator using DetS as first stage estimator.
+
+    Note: The tuning parameter for second stage M estimator is currently only
+    available for a small number of variables and only three values of
+    efficiency. For other cases, the user has to provide the norm instance
+    with desirec tuning parameter.
+
+    Parameters
+    ----------
+    data : array-like
+        Multivariate data set with observation in rows and variables in
+        columns.
+    norm : norm instance
+        If None, then TukeyBiweight norm is used.
+        (Currently no other norms are supported for calling the initial
+        S-estimator)
+        If ``norm`` is an instance of TukeyBiweight, then it will be used in
+        the second stage M-estimation. The ``efficiency`` argument is ignored
+        and the tuning parameter of the user provided instance is not changed.
+    breakdown_point : float in (0, 0.5]
+        Breakdown point for first stage S-estimator.
+    efficiency : float
+        Asymptotic efficiency of second stage M estimator.
+
+    Notes
+    -----
+    Current limitation is that only TukeyBiweight is supported.
+
+    The tuning parameter for second stage M estimator uses a table of values
+    for number of variables up to 15 and efficiency in
+    [0.75, 0.8, 0.85, 0.9, 0.95, 0.975, 0.99]. The tuning parameter for other
+    cases needs to be computed by numerical integration and rootfinding.
+    Alternatively, the user can provide a norm instance with desired tuning
+    parameter.
+
+    References
+    ----------
+    ..[1] Hubert, Mia, Peter Rousseeuw, Dina Vanpaemel, and Tim Verdonck. 2015.
+       “The DetS and DetMM Estimators for Multivariate Location and Scatter.”
+       Computational Statistics & Data Analysis 81 (January): 64–75.
+       https://doi.org/10.1016/j.csda.2014.07.013.
+    ..[2] Hubert, Mia, Peter J. Rousseeuw, and Tim Verdonck. 2012. “A
+       Deterministic Algorithm for Robust Location and Scatter.” Journal of
+       Computational and Graphical Statistics 21 (3): 618–37.
+       https://doi.org/10.1080/10618600.2012.672100.
+    ..[3] Lopuhaä, Hendrik P. 1989. “On the Relation between S-Estimators and
+       M-Estimators of Multivariate Location and Covariance.” The Annals of
+       Statistics 17 (4): 1662–83.
+    ..[4] Salibián-Barrera, Matías, Stefan Van Aelst, and Gert Willems. 2006.
+       “Principal Components Analysis Based on Multivariate MM Estimators with
+       Fast and Robust Bootstrap.” Journal of the American Statistical
+       Association 101 (475): 1198–1211.
+    ..[5] Tatsuoka, Kay S., and David E. Tyler. 2000. “On the Uniqueness of
+       S-Functionals and M-Functionals under Nonelliptical Distributions.” The
+       Annals of Statistics 28 (4): 1219–43.
+    """
+
+    def __init__(self, data, norm=None, breakdown_point=0.5, efficiency=0.95):
+        self.data = np.asarray(data)
+        self.nobs, k_vars = self.data.shape
+        self.k_vars = k_vars
+        self.breakdown_point = breakdown_point
+
+        # self.scale_bias = scale_bias
+        if norm is None:
+            norm = rnorms.TukeyBiweight()
+            c = rtools.tukeybiweight_mvmean_eff(k_vars, efficiency)
+            norm._set_tuning_param(c, inplace=True)
+            # scale_bias is not used for second stage MM norm
+            # self.scale_bias = rtools.scale_bias_cov_biw(c, k_vars)[0]
+        elif not isinstance(norm, rnorms.TukeyBiweight):
+            raise NotImplementedError("only Biweight norm is supported")
+        # We allow tukeybiweight norm instance with user provided c
+
+        self.norm = norm
+        # model for second stage M-estimator
+        self.mod = CovM(data, norm_mean=norm, norm_scatter=norm,
+                        scale_bias=None, method="S")
+
+    def fit(self, maxiter=100):
+        """Estimate model parameters.
+
+        Parameters
+        ----------
+        maxiter : int
+            Maximum number of iterations in the second stage M-estimation.
+        fit args : dict
+            currently missing
+
+        Returns
+        -------
+        Instance of a results or holder class.
+
+        Notes
+        -----
+        This uses CovDetS for the first stage estimation and CovM with fixed
+        scale in the second stage MM-estimation.
+
+        TODO: fit options are missing.
+
+        """
+        # first stage estimate
+        mod_s = CovDetS(
+            self.data,
+            norm=None,
+            breakdown_point=self.breakdown_point
+            )
+        res_s = mod_s.fit()
+
+        res = self.mod.fit(
+            start_mean=res_s.mean,
+            start_shape=res_s.shape,
+            start_scale=res_s.scale,
+            maxiter=maxiter,
+            update_scale=False,
+            )
+
+        return res
diff --git a/statsmodels/robust/norms.py b/statsmodels/robust/norms.py
--- a/statsmodels/robust/norms.py
+++ b/statsmodels/robust/norms.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-# TODO: add plots to weighting functions for online docs.
+from . import tools as rtools
 
 
 def _cabs(x):
@@ -41,6 +41,9 @@ class RobustNorm:
 
     continuous = 1
 
+    def __repr__(self):
+        return self.__class__.__name__
+
     def rho(self, z):
         """
         The robust criterion estimator function.
@@ -199,12 +202,16 @@ class HuberT(RobustNorm):
     def __init__(self, t=1.345):
         self.t = t
 
-    def _set_tuning_param(self, c):
+    def _set_tuning_param(self, c, inplace=False):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        self.t = c
+        if inplace:
+            self.t = c
+            return self
+        else:
+            return self.__class__(t=c)
 
     def max_rho(self):
         return np.inf
@@ -322,6 +329,18 @@ class RamsayE(RobustNorm):
     def __init__(self, a=.3):
         self.a = a
 
+    def _set_tuning_param(self, c, inplace=False):
+        """Set and change the tuning parameter of the Norm.
+
+        Warning: this needs to wipe cached attributes that depend on the param.
+        """
+        # todo : change default to inplace=False, when tools are fixed
+        if inplace:
+            self.a = c
+            return self
+        else:
+            return self.__class__(a=c)
+
     def max_rho(self):
         return np.inf
 
@@ -419,12 +438,16 @@ class AndrewWave(RobustNorm):
     def __init__(self, a=1.339):
         self.a = a
 
-    def _set_tuning_param(self, a):
+    def _set_tuning_param(self, c, inplace=False):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        self.a = a
+        if inplace:
+            self.a = c
+            return self
+        else:
+            return self.__class__(a=c)
 
     def max_rho(self):
         return 2 * self.a**2
@@ -553,12 +576,16 @@ class TrimmedMean(RobustNorm):
     def __init__(self, c=2.):
         self.c = c
 
-    def _set_tuning_param(self, c):
+    def _set_tuning_param(self, c, inplace=False):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        self.c = c
+        if inplace:
+            self.c = c
+            return self
+        else:
+            return self.__class__(c=c)
 
     def max_rho(self):
         return self.rho(self.c)
@@ -674,14 +701,20 @@ def __init__(self, a=2., b=4., c=8.):
         self.b = b
         self.c = c
 
-    def _set_tuning_param(self, c):
+    def _set_tuning_param(self, c, inplace=False):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        self.c = c
-        self.a = c / 4
-        self.b = c / 2
+        a = c / 4
+        b = c / 2
+        if inplace:
+            self.c = c
+            self.a = a
+            self.b = b
+            return self
+        else:
+            return self.__class__(a=a, b=b, c=c)
 
     def max_rho(self):
         return self.rho(self.c)
@@ -862,12 +895,49 @@ class TukeyBiweight(RobustNorm):
     def __init__(self, c=4.685):
         self.c = c
 
-    def _set_tuning_param(self, c):
+    def __repr__(self):
+        return f"{self.__class__.__name__}(c={self.c})"
+
+    @classmethod
+    def get_tuning(cls, bp=None, eff=None):
+        """Tuning parameter for given breakdown point or efficiency.
+
+        This currently only return values from a table.
+
+        Parameters
+        ----------
+        bp : float in [0.05, 0.5] or None
+            Required breakdown point
+            Either bp or eff has to be specified, but not both.
+        eff : float or None
+            Required asymptotic efficiency.
+            Either bp or eff has to be specified, but not both.
+
+        Returns
+        -------
+        float : tuning parameter.
+
+        """
+        if ((bp is None and eff is None) or
+                (bp is not None and eff is not None)):
+            raise ValueError("exactly one of bp and eff needs to be provided")
+
+        if bp is not None:
+            return rtools.tukeybiweight_bp[bp]
+        elif eff is not None:
+            return rtools.tukeybiweight_eff[eff]
+
+    def _set_tuning_param(self, c, inplace=False):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        self.c = c
+        # todo : change default to inplace=False, when tools are fixed
+        if inplace:
+            self.c = c
+            return self
+        else:
+            return self.__class__(c=c)
 
     def max_rho(self):
         return self.rho(self.c)
@@ -891,7 +961,7 @@ def rho(self, z):
         Returns
         -------
         rho : ndarray
-            rho(z) = -(1 - (z/c)**2)**3 * c**2/6.   for \|z\| <= R
+            rho(z) = -(1 - (z/c)**2)**3 * c**2/6 + c**2/6   for \|z\| <= R
 
             rho(z) = 0                              for \|z\| > R
         """
@@ -983,12 +1053,16 @@ def __init__(self, c=3.61752, k=4):
         self.c = c
         self.k = k
 
-    def _set_tuning_param(self, c):
+    def _set_tuning_param(self, c, inplace=False):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        self.c = c
+        if inplace:
+            self.c = c
+            return self
+        else:
+            return self.__class__(c=c, k=self.k)
 
     def max_rho(self):
         return self.rho(self.c)
@@ -1121,12 +1195,16 @@ def __init__(self, c=2.3849, df=4):
         self.c = c
         self.df = df
 
-    def _set_tuning_param(self, c):
+    def _set_tuning_param(self, c, inplace=False):
         """Set and change the tuning parameter of the Norm.
 
         Warning: this needs to wipe cached attributes that depend on the param.
         """
-        self.c = c
+        if inplace:
+            self.c = c
+            return self
+        else:
+            return self.__class__(c=c, df=self.df)
 
     def max_rho(self):
         return np.inf
diff --git a/statsmodels/robust/resistant_linear_model.py b/statsmodels/robust/resistant_linear_model.py
new file mode 100644
--- /dev/null
+++ b/statsmodels/robust/resistant_linear_model.py
@@ -0,0 +1,248 @@
+"""
+Created on Apr. 19, 2024 12:17:03 p.m.
+
+Author: Josef Perktold
+License: BSD-3
+"""
+
+import numpy as np
+from statsmodels.base.model import Model
+from statsmodels.regression.linear_model import OLS
+from statsmodels.tools.testing import Holder
+
+from statsmodels.robust.robust_linear_model import RLM
+import statsmodels.robust.norms as rnorms
+import statsmodels.robust.scale as rscale
+from statsmodels.robust.covariance import _get_detcov_startidx
+
+
+class RLMDetS(Model):
+    """S-estimator for linear model with deterministic starts.
+
+
+    Notes
+    -----
+    This estimator combines the method of Fast-S regression (Saliban-Barrera
+    et al 2006) using starting sets similar to the deterministic estimation of
+    multivariate location and scatter DetS and DetMM of Hubert et al (2012).
+
+
+    References
+    ----------
+
+    .. [1] Hubert, Mia, Peter J. Rousseeuw, and Tim Verdonck. 2012. “A
+       Deterministic Algorithm for Robust Location and Scatter.” Journal of
+       Computational and Graphical Statistics 21 (3): 618–37.
+       https://doi.org/10.1080/10618600.2012.672100.
+
+    .. [2] Hubert, Mia, Peter Rousseeuw, Dina Vanpaemel, and Tim Verdonck.
+       2015. “The DetS and DetMM Estimators for Multivariate Location and
+       Scatter.” Computational Statistics & Data Analysis 81 (January): 64–75.
+       https://doi.org/10.1016/j.csda.2014.07.013.
+
+    .. [3] Rousseeuw, Peter J., Stefan Van Aelst, Katrien Van Driessen, and
+        Jose Agulló. 2004. “Robust Multivariate Regression.”
+        Technometrics 46 (3): 293–305.
+
+    .. [4] Salibian-Barrera, Matías, and Víctor J. Yohai. 2006. “A Fast
+       Algorithm for S-Regression Estimates.” Journal of Computational and
+       Graphical Statistics 15 (2): 414–27.
+
+
+    """
+
+    def __init__(self, endog, exog, norm=None, breakdown_point=0.5,
+                 col_indices=None, include_endog=False):
+        super().__init__(endog, exog)
+
+        if norm is None:
+            norm = rnorms.TukeyBiweight()
+
+        tune = norm.get_tuning(bp=breakdown_point)
+        c = tune[0]
+        scale_bias = tune[2]
+        norm = norm._set_tuning_param(c, inplace=False)
+        self.mscale = rscale.MScale(norm, scale_bias)
+
+        self.norm = norm
+        self.breakdown_point = breakdown_point
+
+        # TODO: detect constant
+        if col_indices is None:
+            exog_start = self.exog[:, 1:]
+        else:
+            exog_start = self.exog[:, col_indices]
+
+        # data for robust mahalanobis distance of starting sets
+        if include_endog:
+            self.data_start = np.column_stack((endog, exog_start))
+        else:
+            self.data_start = exog_start
+
+    def _get_start_params(self, h):
+        # I think we should use iterator with yield
+
+        if self.data_start.shape[1] == 0 and self.exog.shape[1] == 1:
+            quantiles = np.quantile(self.endog, [0.25, 0.5, 0.75])
+            start_params_all = [np.atleast_1d([q]) for q in quantiles]
+            return start_params_all
+
+        starts = _get_detcov_startidx(
+            self.data_start, h, options_start=None, methods_cov="all")
+
+        start_params_all = [
+            OLS(self.endog[idx], self.exog[idx]).fit().params
+            for (idx, method) in starts
+            ]
+        return start_params_all
+
+    def _fit_one(self, start_params, maxiter=100):
+        mod = RLM(self.endog, self.exog, M=self.norm)
+        res = mod.fit(start_params=start_params,
+                      scale_est=self.mscale,
+                      maxiter=maxiter)
+        return res
+
+    def fit(self, h, maxiter=100, maxiter_step=5, start_params_extra=None):
+
+        start_params_all = self._get_start_params(h)
+        if start_params_extra:
+            start_params_all.extend(start_params_extra)
+        res = {}
+        for ii, sp in enumerate(start_params_all):
+            res_ii = self._fit_one(sp, maxiter=maxiter_step)
+            res[ii] = Holder(
+                scale=res_ii.scale,
+                params=res_ii.params,
+                method=ii,  # method  # TODO need start set method
+                )
+
+        scale_all = np.array([i.scale for i in res.values()])
+        scale_sorted = np.argsort(scale_all)
+        best_idx = scale_sorted[0]
+
+        # TODO: iterate until convergence if start fits are not converged
+        res_best = self._fit_one(res[best_idx].params, maxiter=maxiter)
+
+        # TODO: add extra start and convergence info
+        res_best._results.results_iter = res
+        # results instance of _fit_once has RLM as `model`
+        res_best.model_dets = self
+        return res_best
+
+
+class RLMDetSMM(RLMDetS):
+    """MM-estimator with S-estimator starting values.
+
+    Parameters
+    ----------
+    endog : array-like, 1-dim
+        Dependent, endogenous variable.
+    exog array-like, 1-dim
+        Inependent, exogenous regressor variables.
+    norm : robust norm
+        Redescending robust norm used for S- and MM-estimation.
+        Default is TukeyBiweight.
+    efficiency : float in (0, 1)
+        Asymptotic efficiency of the MM-estimator (used in second stage).
+    breakdown_point : float in (0, 0.5)
+        Breakdown point of the preliminary S-estimator.
+    col_indices : None or array-like of ints
+        Index of columns of exog to use in the mahalanobis distance computation
+        for the starting sets of the S-estimator.
+        Default is all exog except first column (constant). Todo: will change
+        when we autodetect the constant column
+    include_endog : bool
+        If true, then the endog variable is combined with the exog variables
+        to compute the the mahalanobis distances for the starting sets of the
+        S-estimator.
+
+    """
+    def __init__(self, endog, exog, norm=None, efficiency=0.95,
+                 breakdown_point=0.5, col_indices=None, include_endog=False):
+        super().__init__(
+            endog,
+            exog,
+            norm=norm,
+            breakdown_point=breakdown_point,
+            col_indices=col_indices,
+            include_endog=include_endog
+            )
+
+        self.efficiency = efficiency
+        if norm is None:
+            norm = rnorms.TukeyBiweight()
+
+        c = norm.get_tuning(eff=efficiency)[0]
+        norm = norm._set_tuning_param(c, inplace=False)
+        self.norm_mean = norm
+
+    def fit(self, h=None, scale_binding=False, start=None):
+        """Estimate the model
+
+        Parameters
+        ----------
+        h : int
+            The size of the initial sets for the S-estimator.
+            Default is ....  (todo)
+        scale_binding : bool
+            If true, then the scale is fixed in the second stage M-estimation,
+            i.e. this is the MM-estimator.
+            If false, then the high breakdown point M-scale is used also in the
+            second stage M-estimation if that estimated scale is smaller than
+            the scale of the preliminary, first stage S-estimato.
+        start : tuple or None
+            If None, then the starting parameters and scale for the second
+            stage M-estimation are taken from the fist stage S-estimator.
+            Alternatively, the starting parameters and starting scale can be
+            provided by the user as tuple (start_params, start_scale). In this
+            case the first stage S-estimation in skipped.
+        maxiter, other optimization parameters are still missing (todo)
+
+        Returns
+        -------
+        results instance
+
+        Notes
+        -----
+        If scale_binding is false, then the estimator is a standard
+        MM-estimator with fixed scale in the second stage M-estimation.
+        If scale_binding is true, then the estimator will try to find an
+        estimate with lower M-scale using the same scale-norm rho as in the
+        first stage S-estimator. If the estimated scale, is not smaller than
+        then the scale estimated in the first stage S-estimator, then the
+        fixed scale MM-estimator is returned.
+
+
+        """
+        norm_m = self.norm_mean
+        if start is None:
+            res_s = super().fit(h)
+            start_params = np.asarray(res_s.params)
+            start_scale = res_s.scale
+        else:
+            start_params, start_scale = start
+            res_s = None
+
+        mod_m = RLM(self.endog, self.exog, M=norm_m)
+        res_mm = mod_m.fit(
+            start_params=start_params,
+            start_scale=start_scale,
+            update_scale=False
+            )
+
+        if not scale_binding:
+            # we can compute this first and skip MM if scale decrease
+            mod_sm = RLM(self.endog, self.exog, M=norm_m)
+            res_sm = mod_sm.fit(
+                start_params=start_params,
+                scale_est=self.mscale
+                )
+
+        if not scale_binding and res_sm.scale < res_mm.scale:
+            res = res_sm
+        else:
+            res = res_mm
+
+        res._results.results_dets = res_s
+        return res
diff --git a/statsmodels/robust/robust_linear_model.py b/statsmodels/robust/robust_linear_model.py
--- a/statsmodels/robust/robust_linear_model.py
+++ b/statsmodels/robust/robust_linear_model.py
@@ -263,6 +263,7 @@ def fit(self, maxiter=50, tol=1e-8, scale_est='mad', init=None, cov='H1',
             wls_results = lm.WLS(self.endog, self.exog).fit()
         else:
             start_params = np.asarray(start_params, dtype=np.double).squeeze()
+            start_params = np.atleast_1d(start_params)
             if (start_params.shape[0] != self.exog.shape[1] or
                     start_params.ndim != 1):
                 raise ValueError('start_params must by a 1-d array with {} '
@@ -276,6 +277,8 @@ def fit(self, maxiter=50, tol=1e-8, scale_est='mad', init=None, cov='H1',
             self.scale = self._estimate_scale(wls_results.resid)
         elif start_scale:
             self.scale = start_scale
+            if not update_scale:
+                self.scale_est = scale_est = "fixed"
 
         history = dict(params=[np.inf], scale=[])
         if conv == 'coefs':
diff --git a/statsmodels/robust/scale.py b/statsmodels/robust/scale.py
--- a/statsmodels/robust/scale.py
+++ b/statsmodels/robust/scale.py
@@ -426,6 +426,9 @@ def __init__(self, chi_func, scale_bias):
         self.chi_func = chi_func
         self.scale_bias = scale_bias
 
+    def __repr__(self):
+        return repr(self.chi_func)
+
     def __call__(self, data, **kwds):
         return self.fit(data, **kwds)
 
diff --git a/statsmodels/robust/tools.py b/statsmodels/robust/tools.py
--- a/statsmodels/robust/tools.py
+++ b/statsmodels/robust/tools.py
@@ -4,7 +4,7 @@
 Author: Josef Perktold
 License: BSD-3
 """
-
+# flake8: noqa E731
 import numpy as np
 from scipy import stats, integrate, optimize
 
@@ -57,8 +57,8 @@ def _var_normal(norm):
 
 
     """
-    num = stats.norm.expect(lambda x: norm.psi(x)**2)
-    denom = stats.norm.expect(lambda x: norm.psi_deriv(x))**2
+    num = stats.norm.expect(lambda x: norm.psi(x)**2)  #noqa E731
+    denom = stats.norm.expect(lambda x: norm.psi_deriv(x))**2  #noqa E731
     return num / denom
 
 
@@ -100,7 +100,7 @@ def _var_normal_jump(norm):
 
 
     """
-    num = stats.norm.expect(lambda x: norm.psi(x)**2)
+    num = stats.norm.expect(lambda x: norm.psi(x)**2)  #noqa E731
 
     def func(x):
         # derivative normal pdf
@@ -142,7 +142,6 @@ def _get_tuning_param(norm, eff, kwd="c", kwargs=None, use_jump=False,
         efficiency.
 
     """
-    kwds = {} if kwargs is None else kwargs
     if bracket is None:
         bracket = [0.1, 10]
 
@@ -150,12 +149,12 @@ def _get_tuning_param(norm, eff, kwd="c", kwargs=None, use_jump=False,
         def func(c):
             # kwds.update({kwd: c})
             # return _var_normal(norm(**kwds)) - 1 / eff
-            norm._set_tuning_param(c)
+            norm._set_tuning_param(c, inplace=True)
             return _var_normal(norm) - 1 / eff
     else:
         def func(c):
-            norm._set_tuning_param(c)
-            return _var_normal_jump(norm(**kwds) - 1 / eff)
+            norm._set_tuning_param(c, inplace=True)
+            return _var_normal_jump(norm) - 1 / eff
 
     res = optimize.brentq(func, *bracket)
     return res
@@ -215,14 +214,15 @@ def tuning_s_estimator_mean(norm, breakdown=None):
 
     def func(c):
         norm_ = norm
-        norm_._set_tuning_param(c)
-        bp = stats.norm.expect(lambda x : norm_.rho(x)) / norm_.max_rho()
+        norm_._set_tuning_param(c, inplace=True)
+        bp = (stats.norm.expect(lambda x : norm_.rho(x)) /  #noqa E731
+              norm_.max_rho())
         return bp
 
     res = []
     for bp in bps:
-        c_bp = optimize.brentq(lambda c0: func(c0) - bp, 0.1, 10)
-        norm._set_tuning_param(c_bp)  # inplace modification
+        c_bp = optimize.brentq(lambda c0: func(c0) - bp, 0.1, 10)  #noqa E731
+        norm._set_tuning_param(c_bp, inplace=True)  # inplace modification
         eff = 1 / _var_normal(norm)
         b = stats.norm.expect(lambda x : norm.rho(x))
         res.append([bp, eff, c_bp, b])
@@ -242,3 +242,321 @@ def func(c):
         )
 
     return res2
+
+
+def scale_bias_cov_biw(c, k_vars):
+    """Multivariate scale bias correction for TukeyBiweight norm.
+
+    This uses the chisquare distribution as reference distribution for the
+    squared Mahalanobis distance.
+    """
+    p = k_vars  # alias for formula
+    chip, chip2, chip4, chip6 = stats.chi2.cdf(c**2, [p, p + 2, p + 4, p + 6])
+    b = p / 2 * chip2 - p * (p + 2) / (2 * c**2) * chip4
+    b += p * (p + 2) * (p + 4) / (6 * c**4) * chip6 + c**2 / 6 * (1 - chip)
+    return b, b / (c**2 / 6)
+
+
+def scale_bias_cov(norm, k_vars):
+    """Multivariate scale bias correction.
+
+
+    Parameter
+    ---------
+    norm : norm instance
+        The rho function of the norm is used in the moment condition for
+        estimating scale.
+    k_vars : int
+        Number of random variables in the multivariate data.
+
+    Returns
+    -------
+    scale_bias: float
+    breakdown_point : float
+        Breakdown point computed as scale bias divided by max rho.
+    """
+
+    rho = lambda x: (norm.rho(np.sqrt(x)))  # noqa
+    scale_bias = stats.chi2.expect(rho, args=(k_vars,))
+    return scale_bias, scale_bias / norm.max_rho()
+
+
+def tuning_s_cov(norm, k_vars, breakdown_point=0.5, limits=()):
+    """Tuning parameter for multivariate S-estimator given breakdown point.
+    """
+    from .norms import TukeyBiweight  # avoid circular import
+
+    if not limits:
+        limits = (0.5, 30)
+
+    if isinstance(norm, TukeyBiweight):
+        def func(c):
+            return scale_bias_cov_biw(c, k_vars)[1] - breakdown_point
+    else:
+        norm = norm._set_tuning_param(2., inplace=False)  # create copy
+
+        def func(c):
+            norm._set_tuning_param(c, inplace=True)
+            return scale_bias_cov(norm, k_vars)[1] - breakdown_point
+
+    p_tune = optimize.brentq(func, limits[0], limits[1])
+    return p_tune
+
+
+def eff_mvmean(norm, k_vars):
+    """Efficiency for M-estimator of multivariate mean at normal distribution.
+
+    This also applies to estimators that are locally equivalent to an
+    M-estimator such as S- and MM-estimators.
+
+    Parameters
+    ----------
+    norm : instance of norm class
+    k_vars : int
+        Number of variables in multivariate random variable, i.e. dimension.
+
+    Returns
+    -------
+    eff : float
+        Asymptotic relative efficiency of mean at normal distribution.
+    alpha : float
+        Numerical integral. Efficiency is beta**2 / alpha
+    beta : float
+        Numerical integral.
+
+    Notes
+    -----
+    This implements equ. (5.3) p. 1671 in Lopuhaä 1989
+
+    References
+    ----------
+
+    .. [1] Lopuhaä, Hendrik P. 1989. “On the Relation between S-Estimators
+       and M-Estimators of Multivariate Location and Covariance.”
+       The Annals of Statistics 17 (4): 1662–83.
+
+    """
+    k = k_vars  # shortcut
+    f_alpha = lambda d: norm.psi(d)**2 / k
+    f_beta = lambda d: ((1 - 1 / k) * norm.weights(d) +
+                        1 / k * norm.psi_deriv(d))
+    alpha = stats.chi(k).expect(f_alpha)
+    beta = stats.chi(k).expect(f_beta)
+    return beta**2 / alpha, alpha, beta
+
+
+def eff_mvshape(norm, k_vars):
+    """Efficiency of M-estimator of multivariate shape at normal distribution.
+
+    This also applies to estimators that are locally equivalent to an
+    M-estimator such as S- and MM-estimators.
+
+    Parameters
+    ----------
+    norm : instance of norm class
+    k_vars : int
+        Number of variables in multivariate random variable, i.e. dimension.
+
+    Returns
+    -------
+    eff : float
+        Asymptotic relative efficiency of mean at normal distribution.
+    alpha : float
+        Numerical integral. Efficiency is beta**2 / alpha
+    beta : float
+        Numerical integral.
+
+    Notes
+    -----
+    This implements sigma_1 in equ. (5.5) p. 1671 in Lopuhaä 1989.
+    Efficiency of shape is approximately 1 / sigma1.
+
+    References
+    ----------
+
+    .. [1] Lopuhaä, Hendrik P. 1989. “On the Relation between S-Estimators
+       and M-Estimators of Multivariate Location and Covariance.”
+       The Annals of Statistics 17 (4): 1662–83.
+
+    """
+
+    k = k_vars  # shortcut
+    f_a = lambda d: k * (k + 2) * norm.psi(d)**2 * d**2
+    f_b = lambda d: norm.psi_deriv(d) * d**2 + (k + 1) * norm.psi(d) * d
+    a = stats.chi(k).expect(f_a)
+    b = stats.chi(k).expect(f_b)
+    return b**2 / a, a, b
+
+
+def tuning_m_cov_eff(norm, k_vars, efficiency=0.95, eff_mean=True, limits=()):
+    """Tuning parameter for multivariate M-estimator given efficiency.
+
+    This also applies to estimators that are locally equivalent to an
+    M-estimator such as S- and MM-estimators.
+
+    Parameters
+    ----------
+    norm : instance of norm class
+    k_vars : int
+        Number of variables in multivariate random variable, i.e. dimension.
+    efficiency : float < 1
+        Desired asymptotic relative efficiency of mean estimator.
+        Default is 0.95.
+    eff_mean : bool
+        If eff_mean is true (default), then tuning parameter is to achieve
+        efficiency of mean estimate.
+        If eff_mean is fale, then tuning parameter is to achieve efficiency
+        of shape estimate.
+    limits : tuple
+        Limits for rootfinding with scipy.optimize.brentq.
+        In some cases the interval limits for rootfinding can be too small
+        and not cover the root. Current default limits are (0.5, 30).
+
+    Returns
+    -------
+    float : Tuning parameter for the norm to achieve desired efficiency.
+        Asymptotic relative efficiency of mean at normal distribution.
+
+    Notes
+    -----
+    This uses numerical integration and rootfinding and will be
+    relatively slow.
+    """
+    if not limits:
+        limits = (0.5, 30)
+
+    # make copy of norm
+    norm = norm._set_tuning_param(1, inplace=False)
+
+    if eff_mean:
+        def func(c):
+            norm._set_tuning_param(c, inplace=True)
+            return eff_mvmean(norm, k_vars)[0] - efficiency
+    else:
+        def func(c):
+            norm._set_tuning_param(c, inplace=True)
+            return eff_mvshape(norm, k_vars)[0] - efficiency
+
+    p_tune = optimize.brentq(func, limits[0], limits[1])
+    return p_tune
+
+
+#  ##### tables
+
+tukeybiweight_bp = {
+    # breakdown point : (tuning parameter, efficiency, scale bias)
+    0.50: (1.547645, 0.286826, 0.199600),
+    0.45: (1.756059, 0.369761, 0.231281),
+    0.40: (1.987965, 0.461886, 0.263467),
+    0.35: (2.251831, 0.560447, 0.295793),
+    0.30: (2.560843, 0.661350, 0.327896),
+    0.25: (2.937015, 0.759040, 0.359419),
+    0.20: (3.420681, 0.846734, 0.390035),
+    0.15: (4.096255, 0.917435, 0.419483),
+    0.10: (5.182361, 0.966162, 0.447614),
+    0.05: (7.545252, 0.992424, 0.474424),
+    }
+
+tukeybiweight_eff = {
+    # efficiency : (tuning parameter, breakdown point)
+    0.65: (2.523102, 0.305646),
+    0.70: (2.697221, 0.280593),
+    0.75: (2.897166, 0.254790),
+    0.80: (3.136909, 0.227597),
+    0.85: (3.443690, 0.197957),
+    0.90: (3.882662, 0.163779),
+    0.95: (4.685065, 0.119414),
+    0.97: (5.596823, 0.087088),
+    0.98: (5.920719, 0.078604),
+    0.99: (7.041392, 0.056969),
+    }
+
+# relative efficiency for M and MM estimator of multivariate location
+# Table 2 from Kudraszow and Maronna JMA 2011
+# k in [1, 2, 3, 4, 5, 10], eff in [0.8, 0.9, 0.95
+# TODO: need to replace with more larger table and more digits.
+# (4, 0.95): 5.76 -. 5.81 to better match R rrcov and numerical integration
+tukeybiweight_mvmean_eff_km = {
+        (1, 0.8): 3.14, (2, 0.8): 3.51, (3, 0.8): 3.82, (4, 0.8): 4.1,
+        (5, 0.8): 4.34, (10, 0.8): 5.39,
+        (1, 0.9): 3.88, (2, 0.9): 4.28, (3, 0.9): 4.62, (4, 0.9): 4.91,
+        (5, 0.9): 5.18, (10, 0.9): 6.38,
+        (1, 0.95): 4.68, (2, 0.95): 5.12, (3, 0.95): 5.48, (4, 0.95): 5.81,
+        (5, 0.95): 6.1, (10, 0.95): 7.67,
+        }
+
+_table_biweight_mvmean_eff = np.array([
+    [2.89717, 3.13691, 3.44369, 3.88266, 4.68506, 5.59682, 7.04139],
+    [3.26396, 3.51006, 3.82643, 4.2821, 5.12299, 6.0869, 7.62344],
+    [3.5721, 3.82354, 4.14794, 4.61754, 5.49025, 6.49697, 8.10889],
+    [3.84155, 4.09757, 4.42889, 4.91044, 5.81032, 6.85346, 8.52956],
+    [4.08323, 4.34327, 4.68065, 5.17267, 6.09627, 7.17117, 8.90335],
+    [4.30385, 4.56746, 4.91023, 5.41157, 6.35622, 7.45933, 9.24141],
+    [4.50783, 4.77466, 5.12228, 5.63199, 6.59558, 7.72408, 9.5512],
+    [4.69828, 4.96802, 5.32006, 5.83738, 6.81817, 7.96977, 9.83797],
+    [4.87744, 5.14986, 5.506, 6.03022, 7.02679, 8.19958, 10.10558],
+    [5.04704, 5.32191, 5.68171, 6.21243, 7.22354, 8.41594, 10.35696],
+    [5.20839, 5.48554, 5.84879, 6.38547, 7.41008, 8.62071, 10.59439],
+    [5.36254, 5.64181, 6.00827, 6.55051, 7.58772, 8.81538, 10.81968],
+    [5.51034, 5.7916, 6.16106, 6.7085, 7.75752, 9.00118, 11.03428],
+    [5.65249, 5.9356, 6.3079, 6.86021, 7.92034, 9.17908, 11.23939],
+    ])
+
+_table_biweight_mvshape_eff = np.array([
+    [3.57210, 3.82354, 4.14794, 4.61754, 5.49025, 6.49697, 8.10889],
+    [3.84155, 4.09757, 4.42889, 4.91044, 5.81032, 6.85346, 8.52956],
+    [4.08323, 4.34327, 4.68065, 5.17267, 6.09627, 7.17117, 8.90335],
+    [4.30385, 4.56746, 4.91023, 5.41157, 6.35622, 7.45933, 9.24141],
+    [4.50783, 4.77466, 5.12228, 5.63199, 6.59558, 7.72408, 9.55120],
+    [4.69828, 4.96802, 5.32006, 5.83738, 6.81817, 7.96977, 9.83797],
+    [4.87744, 5.14986, 5.50600, 6.03022, 7.02679, 8.19958, 10.10558],
+    [5.04704, 5.32191, 5.68171, 6.21243, 7.22354, 8.41594, 10.35696],
+    [5.20839, 5.48554, 5.84879, 6.38547, 7.41008, 8.62071, 10.59439],
+    [5.36254, 5.64181, 6.00827, 6.55051, 7.58772, 8.81538, 10.81968],
+    [5.51034, 5.79160, 6.16106, 6.70849, 7.75752, 9.00118, 11.03428],
+    [5.65249, 5.93560, 6.30790, 6.86021, 7.92034, 9.17908, 11.23939],
+    [5.78957, 6.07443, 6.44938, 7.00630, 8.07692, 9.34991, 11.43603],
+    [5.92206, 6.20858, 6.58604, 7.14731, 8.22785, 9.51437, 11.62502],
+    ])
+
+
+def _convert_to_dict_mvmean_effs(eff_mean=True):
+    effs_mvmean = [0.75, 0.8, 0.85, 0.9, 0.95, 0.975, 0.99]
+    ks = list(range(1, 15))
+    if eff_mean:
+        table = _table_biweight_mvmean_eff
+    else:
+        table = _table_biweight_mvshape_eff
+    tp = {}
+    for i, k in enumerate(ks):
+        for j, eff in enumerate(effs_mvmean):
+            tp[(k, eff)] = table[i, j]
+
+    return tp
+
+
+tukeybiweight_mvmean_eff_d = _convert_to_dict_mvmean_effs(eff_mean=True)
+tukeybiweight_mvshape_eff_d = _convert_to_dict_mvmean_effs(eff_mean=False)
+
+
+def tukeybiweight_mvmean_eff(k, eff, eff_mean=True):
+    """tuning parameter for biweight norm to achieve efficiency for mv-mean.
+
+    Uses values from precomputed table if available, otherwise computes it
+    numerically and adds it to the module global dict.
+
+    """
+    if eff_mean:
+        table_dict = tukeybiweight_mvmean_eff_d
+    else:
+        table_dict = tukeybiweight_mvshape_eff_d
+
+    try:
+        tp = table_dict[(k, eff)]
+    except KeyError:
+        # compute and cache
+        from .norms import TukeyBiweight  # avoid circular import
+        norm = TukeyBiweight(c=1)
+        tp = tuning_m_cov_eff(norm, k, efficiency=eff, eff_mean=eff_mean)
+        table_dict[(k, eff)] = tp
+    return tp
diff --git a/statsmodels/stats/covariance.py b/statsmodels/stats/covariance.py
--- a/statsmodels/stats/covariance.py
+++ b/statsmodels/stats/covariance.py
@@ -149,6 +149,18 @@ def transform_corr_normal(corr, method, return_var=False, possdef=True):
         return corr_n
 
 
+def corr_rank(data):
+    """Spearman rank correlation
+
+    simplified version of scipy.stats.spearmanr
+    """
+    x = np.asarray(data)
+    axisout = 0
+    ar = np.apply_along_axis(stats.rankdata, axisout, x)
+    corr = np.corrcoef(ar, rowvar=False)
+    return corr
+
+
 def corr_normal_scores(data):
     """Gaussian rank (normal scores) correlation
 


Test patch:
diff --git a/statsmodels/robust/tests/test_covariance.py b/statsmodels/robust/tests/test_covariance.py
--- a/statsmodels/robust/tests/test_covariance.py
+++ b/statsmodels/robust/tests/test_covariance.py
@@ -6,6 +6,7 @@
 import pandas as pd
 
 from statsmodels import robust
+import statsmodels.robust.norms as robnorms
 import statsmodels.robust.covariance as robcov
 import statsmodels.robust.scale as robscale
 
@@ -177,6 +178,121 @@ def test_tyler():
     assert res1.n_iter == 55
 
 
+def test_cov_ms():
+    # use CovM as local CovS
+
+    # > CovSest(x)  using package rrcov
+    # same result with > CovSest(x, method="sdet")
+    # but scale difers with method="biweight"
+    mean_r = np.array([1.53420879676, 1.82865741024, 1.65565146981])
+    cov_r = np.array([
+        [1.8090846049573, 0.0479283121828, 0.2446369025717],
+        [0.0479283121828, 1.8189886310494, 0.2513025527579],
+        [0.2446369025717, 0.2513025527579, 1.7287983150484],
+        ])
+
+    scale2_r = np.linalg.det(cov_r) ** (1/3)
+    shape_r = cov_r / scale2_r
+    scale_r = np.sqrt(scale2_r)
+
+    exog_df = dta_hbk[["X1", "X2", "X3"]]
+    mod = robcov.CovM(exog_df)
+    # with default start, default start could get wrong local optimum
+    res = mod.fit()
+    assert_allclose(res.mean, mean_r, rtol=1e-5)
+    assert_allclose(res.shape, shape_r, rtol=1e-5)
+    assert_allclose(res.cov, cov_r, rtol=1e-5)
+    assert_allclose(res.scale, scale_r, rtol=1e-5)
+
+    # with results start
+    res = mod.fit(start_mean=mean_r, start_shape=shape_r, start_scale=scale_r)
+    assert_allclose(res.mean, mean_r, rtol=1e-5)
+    assert_allclose(res.shape, shape_r, rtol=1e-5)
+    assert_allclose(res.cov, cov_r, rtol=1e-5)
+    assert_allclose(res.scale, scale_r, rtol=1e-5)
+
+    mod_s = robcov.CovDetS(exog_df)
+    res = mod_s.fit()
+    assert_allclose(res.mean, mean_r, rtol=1e-5)
+    assert_allclose(res.shape, shape_r, rtol=1e-5)
+    assert_allclose(res.cov, cov_r, rtol=1e-5)
+    assert_allclose(res.scale, scale_r, rtol=1e-5)
+
+
+def test_covdetmcd():
+
+    # results from rrcov
+    # > cdet = CovMcd(x = hbk, raw.only = TRUE, nsamp = "deterministic",
+    #                 use.correction=FALSE)
+    cov_dmcd_r = np.array("""
+    2.2059619213639   0.0223939863695   0.7898958050933   0.4060613360808
+    0.0223939863695   1.1384166802155   0.4315534571891  -0.2344041030201
+    0.7898958050933   0.4315534571891   1.8930117467493  -0.3292893001459
+    0.4060613360808  -0.2344041030201  -0.3292893001459   0.6179686100845
+    """.split(), float).reshape(4, 4)
+
+    mean_dmcd_r = np.array([1.7725, 2.2050, 1.5375, -0.0575])
+
+    mod = robcov.CovDetMCD(dta_hbk)
+    res = mod.fit(40, maxiter_step=100, reweight=False)
+    assert_allclose(res.mean, mean_dmcd_r, rtol=1e-5)
+    assert_allclose(res.cov, cov_dmcd_r, rtol=1e-5)
+
+    # with reweighting
+    # covMcd(x = hbk, nsamp = "deterministic", use.correction = FALSE)
+    # iBest: 5; C-step iterations: 7, 7, 7, 4, 6, 6
+    # Log(Det.):  -2.42931967153
+
+    mean_dmcdw_r = np.array([1.5338983050847, 1.8322033898305, 1.6745762711864,
+                            -0.0728813559322])
+    cov_dmcdw_r = np.array("""
+    1.5677744869295   0.09285770205078   0.252076010128   0.13873444408300
+    0.0928577020508   1.56769177397171   0.224929617385  -0.00516128856542
+    0.2520760101278   0.22492961738467   1.483829106079  -0.20275013775619
+    0.1387344440830  -0.00516128856542  -0.202750137756   0.43326701543885
+    """.split(), float).reshape(4, 4)
+
+    mod = robcov.CovDetMCD(dta_hbk)
+    res = mod.fit(40, maxiter_step=100)  # default is reweight=True
+    assert_allclose(res.mean, mean_dmcdw_r, rtol=1e-5)
+    # R uses different trimming correction
+    # compare only shape (using trace for simplicity)
+    shape = res.cov / np.trace(res.cov)
+    shape_r = cov_dmcdw_r / np.trace(cov_dmcdw_r)
+    assert_allclose(shape, shape_r, rtol=1e-5)
+
+
+def test_covdetmm():
+
+    # results from rrcov
+    # CovMMest(x = hbk, eff.shape=FALSE,
+    #          control=CovControlMMest(sest=CovControlSest(method="sdet")))
+    cov_dmm_r = np.array("""
+        1.72174266670826 0.06925842715939 0.20781848922667 0.10749343153015
+        0.06925842715939 1.74566218886362 0.22161135221404 -0.00517023660647
+        0.20781848922667 0.22161135221404 1.63937749762534 -0.17217102475913
+        0.10749343153015 -0.00517023660647 -0.17217102475913 0.48174480967136
+        """.split(), float).reshape(4, 4)
+
+    mean_dmm_r = np.array([1.5388643420460, 1.8027582110408, 1.6811517253521,
+                           -0.0755069488908])
+
+    # using same c as rrcov
+    c = 5.81031555752526
+    mod = robcov.CovDetMM(dta_hbk, norm=robnorms.TukeyBiweight(c=c))
+    res = mod.fit()
+
+    assert_allclose(res.mean, mean_dmm_r, rtol=1e-5)
+    assert_allclose(res.cov, cov_dmm_r, rtol=1e-5, atol=1e-5)
+
+    # using c from table,
+    mod = robcov.CovDetMM(dta_hbk)
+    res = mod.fit()
+
+    assert_allclose(res.mean, mean_dmm_r, rtol=1e-3)
+    assert_allclose(res.cov, cov_dmm_r, rtol=1e-3, atol=1e-3)
+
+
 def test_robcov_SMOKE():
     # currently only smoke test or very loose comparisons to dgp
     nobs, k_vars = 100, 3
diff --git a/statsmodels/robust/tests/test_tools.py b/statsmodels/robust/tests/test_tools.py
--- a/statsmodels/robust/tests/test_tools.py
+++ b/statsmodels/robust/tests/test_tools.py
@@ -49,24 +49,26 @@
 @pytest.mark.parametrize("case", results_menenez)
 def test_eff(case):
     norm, res2 = case
+    use_jump = False
 
     if norm.continuous == 2:
         var_func = _var_normal
     else:
         var_func = _var_normal_jump
+        use_jump = True
 
     res_eff = []
     for c in res2:
-        norm._set_tuning_param(c)
+        norm._set_tuning_param(c, inplace=True)
         res_eff.append(1 / var_func(norm))
 
     assert_allclose(res_eff, effs, atol=0.0005)
 
     for c in res2:
         # bp = stats.norm.expect(lambda x : norm.rho(x)) / norm.rho(norm.c)
-        norm._set_tuning_param(c)
-        eff = 1 / _var_normal(norm)
-        tune = _get_tuning_param(norm, eff)
+        norm._set_tuning_param(c, inplace=True)
+        eff = 1 / var_func(norm)
+        tune = _get_tuning_param(norm, eff, use_jump=use_jump)
         assert_allclose(tune, c, rtol=1e-6, atol=5e-4)