determine constant for biweight (bisquare) function

Question

In robust statistics a biweight (bisquare) function is defined as follows

$$\rho \left( x \right) = \gamma\left( {1 - {{\left( {1 - {{\left( {\frac{x}{c}} \right)}^2}} \right)}^3}} \right){{\bf{1}}_{\left| x \right| \le c}} + \gamma{{\bf{1}}_{\left| x \right| > c}}$$

For an n-dimensional random variable the constant $\gamma$ can be determined using the consistency equation and c is given by the quantile of ${\chi}^2_{1-\alpha,n}$ $$ E\left[ {\rho \left( \|X \|_2\right)} \right] = n$$

If $X$ is standard n-dimensional multivariate random variable then we need to solve the following to get the value of $\gamma$ $$\int ... \int_{\textbf{R}^n} {\rho \left( {\sqrt {{z^ \top }z} } \right)} {\frac{1}{{{{\left( {2\pi } \right)}^{n/2}}}}\exp \left\{ { - \frac{1}{2}{z^ \top }z} \right\}}\,dz_n\,dz_{n-1}\,...dz_1 = n$$

This integral can be solved numerically for low dimensions or by simulation by simulating variables from a standard multivariate normal distribution, calculating the average loss using the bisquare function and setting up an optimization problem to minimize the absolute difference between the average loss and dimension n.

However, I am also interested in evaluating the integral analytically or atleast make the integral more tractable to solve. The idea I have is that, the sum of squares of normal random variables has a chi-square distribution. Hence a transformation such as, $k = {z^ \top }z$, where k is chi-square random variable can collapse the integral to a univariate integral. Can someone point an example on how to convert a multidimensional integral to univariate integral or point out what the integrator and limits of the new integrator in the transformed space should be.

Answer 1

I will essentially summarize what is said in [1,p185-186]. Your notation is in many places different from the usual one, so I took the liberty of rewriting some of the introductory concepts so that we're on the same tune.

Denote $w(d(\pmb x))=\rho^\prime(d(\pmb x))$, the weight function corresponding to the (multivariate version of the) bi-square loss function:

$$\rho^\prime(d(\pmb x))=1_{|d(\pmb x)|\leq c}\left(d(\pmb x)\left(1-\left(\frac{d(\pmb x)}{c}\right)^2\right)^2\right)$$

and $d(\pmb x)=(\pmb x-\hat{\pmb u}_{\infty})^\top\hat{\pmb\varSigma}_{\infty}^{-1}(\pmb x-\hat{\pmb u}_{\infty})$. Then we want the final estimates to be consistent in the sense that:

$$(0)\quad E\left(w(d(\pmb x))(\pmb x-\hat{\pmb u}_{\infty})(\pmb x-\hat{\pmb u}_{\infty})^\top\right)=\hat{\pmb\varSigma}_{\infty}$$

where $\pmb\mu$ is the parameter, $\hat{\pmb\mu}$ is the estimate of $\pmb\mu$ and $\hat{\pmb\mu}_{\infty}$ the asymptotic value of said estimate (and likewise for $\hat{\pmb\varSigma}_{\infty}$).

The multivariate S estimator is affine equivariant, so we can set w.l.o.g. set $\hat{\pmb\varSigma}_{\infty}=\pmb I_n$ (the rank $n$ identity matrix) and $\hat{\pmb u}_{\infty}=\pmb 0_n$ (the $n$ dimensional origin). So $d(\pmb x)$ becomes $||\pmb x||^2/c$ and $(0)$ simplifies to:

$$(1)\quad E\left(w\left(\frac{||\pmb x||^2}{c}\right)||\pmb x||^2\right)=nc.$$

Since as you pointed out, when $\pmb x$ is $n$-variate standard normal $||\pmb x||^2$ has a $\chi^2_n$ distribution, we obtain a consistent estimate of $\pmb\varSigma$ in the normal case by replacing $\hat{\pmb\varSigma}_n$ by $\hat{\pmb\varSigma}_n/c$ with $c$ defined as the solution of:

$$(2)\quad\int_0^{\infty}w\left(\frac{z}{c}\right)\frac{z}{c}g(z)dz=n$$

where $g$ is the density of the $\chi_n^2$ distribution. I guess you could use a root finding algorithm and numerical integration to solve $(2)$ for $c$, depending on the application it could be fairly quick nowadays.

---

Edit:

whuber's comment is spot on in that the problem above can be simplified considerably using the first 6 moments of the $\chi^2_n$ distribution (see [2]). It can also be done using the first 3 moments of a gamma distribution, but it turns out this second solution is more expensive computationally (see [2] for the source of this claim). The details are well explained in [2], but for completeness, I will put the main result. The constant $c$ in $(2)$ must satisfy [2,equation 32]:

$$(3)\quad\frac{1}{2}=\frac{3n}{c^2}\left(\Psi_{n+2}^c-(n+2)\frac{\Psi_{n+4}^c}{c^2}+(n+2)(n+4)\frac{\Psi_{n+6}^c}{3c^4}+\frac{c^2}{3}\frac{1-\Psi_{n}^c}{n}\right)$$

where $\Psi_n^c=\text{Pr}(\chi^2_n<c^2)$ and the $\frac{1}{2}$ on the lhs of $(3)$ is actually the breakdown point you want your estimate to be consist at (i.e. a number in $(0,0.5]$). So compared to $(2)$, $(3)$ bypasses the need for numerical integration so that it can be solved using a root finding algorithm only, which is of course much faster.

---

• [1] Robust Statistics: Theory and Methods (2006). Maronna, R. A., Martin, D. R. and V. J. Yohai.
• [2] On consistency factors and efficiency of robust S-estimators (2014). M. Riani , A. Cerioli and F. Torti.