3
$\begingroup$

Was wondering if anyone knows of an R package to estimate the Cauchy-M estimator of regression (see for example the end of this section, but with simultaneous estimation of the scale parameter as in section 2 of (1)).

(1) Mizera, I. Müller, C. H. (2002). Breakdown points of Cauchy regression-scale estimators, Statistics & Probability Letters, Volume 57, Issue 1, Pages 79-89.

$\endgroup$
5
  • $\begingroup$ MASS::rlm will let you specify a $\psi$ function, though I think the choice of scale estimate is limited to a few options. $\endgroup$
    – Glen_b
    Apr 29 '13 at 11:23
  • $\begingroup$ Outside of that, there are some functions in package robustbase that might let you get near to what you want, but I'm not certain of that. $\endgroup$
    – Glen_b
    Apr 29 '13 at 11:41
  • $\begingroup$ @Glen_b: thanks. Actually, I believe this problem is somewhat more popular in the context of Bayesian estimation. I was hoping that someone more familiar with Bayesian estimation could help but I notice smillig disagrees with that assessment. $\endgroup$
    – user603
    Apr 29 '13 at 11:54
  • 1
    $\begingroup$ A very straight-forward approach is via the EM algorithm, as Cauchy-likelihood is a scale mixture of normal likelihoods. So you get a weighted least squares regression were the weights are adaptively chosen based on the "current" residual. $\endgroup$ Apr 29 '13 at 13:33
  • $\begingroup$ @probabilityislogic: I think I understand you, but can you suggest a practical way to hack a EM-algorithm for regression to achieve this? $\endgroup$
    – user603
    Apr 29 '13 at 14:00
5
$\begingroup$

This arxiv paper may be of use to you:

Data augmentation for non-Gaussian regression models using variance-mean mixtures

Very simple to code up and implement, you could use the standard linear model function lm(y~formula,weights=w) for the "M" step (including estimation of $\sigma$) where the weights are taken from the "E" step.

Note that for a Cauchy RV we have that $(e_i|\sigma^2\omega_i)\sim N(0,\sigma^2\omega_i^{-1})$ and $\omega_i\sim \chi^2(1)$ implies $(e_i|\sigma^2)\sim Cauchy(0,\sigma^2)$. The derivative of the negative log likelihood is $\frac{2e_i}{\sigma^2+e_i^2}$ which is $f'()$ in proposition of the paper. This means that the "E" step is given as:

$$\hat{\omega}_i=\frac{2\sigma^2}{\sigma^2+\hat{r}_i^2}$$

where $\hat{r}_i=y_i-x_i^T\hat{\beta}$ is the current residual (easily obtained from lm()$residuals in the previous "M" step).

The tricky part with this is to note that the second derivative of the negative log likelihood is given as $$\frac{2}{\sigma^2+e_i^2}-\frac{4e_i^2}{(\sigma^2+e_i^2)^2}$$ which mean the optimisation problem is not convex - so you may have multiple modes. EM will surely find one of them, but you need to do re-starts at different values to make sure you have the actual maximum.

UPDATE

As mentioned in my comment, I have added some R code demonstrating the multi-modal nature of Cauchy regression. This is especially true in just the cases where Cauchy regression is most useful - when there are outliers that influence the regression fit. This example shows this by using a point of high leverage which doesn't follow the rest of the data.

#em algorithm for Cauchy regression
cauchylm <- function(y,x,beta_init=NULL){
 n <- length(y)
 if(is.null(beta_init)){
 ols_reg <- lm(y~x)
 beta_old <- ols_reg$coefficients
     sigma_old <- summary(ols_reg)$sigma
 }else{
 beta_old <- beta_init
 ols_reg <- list(residuals=y-beta_init[1]-x*beta_init[2])
 sigma_old <- sqrt(sum(ols_reg$residuals^2)/(n-2))
     }
     iter<- 1
     prec <- 1
     while((iter<1000)&(prec>0.00000001)){
     weights <- 2*sigma_old^2/(sigma_old^2+ols_reg$residuals^2)
     ols_reg <- lm(y~x,weights=weights)
     beta <- ols_reg$coefficients
     sigma <- summary(ols_reg)$sigma
 prec <- max(abs(c(beta,sigma)-c(beta_old,sigma_old))/(abs(c(beta_old,sigma_old))+0.001))
 beta_old <- beta
 sigma_old <- sigma
 iter <- iter+1
 }
return(list(beta=beta,sigma=sigma,ols_reg=ols_reg,iter=iter,prec=prec))
}

#this example creates a point of high leverage, showing the multimodality

x<- runif(10)
beta <- -1
alpha <- 1
y <- alpha + beta*x + 0.3*rnorm(10)

x0 <- 4
y0 <- x0+0.3*rnorm(1)

# Cauchy regression using OLS estimates as starting values
reg1 <- cauchylm(c(y,y0),c(x,x0))
# Cauchy regression using true values as starting values
reg2 <- cauchylm(c(y,y0),c(x,x0),c(alpha,beta))
# Cauchy regression using random values as starting values
reg3 <- cauchylm(c(y,y0),c(x,x0),rnorm(2))
# Cauchy regression using random values as starting values
reg4 <- cauchylm(c(y,y0),c(x,x0),c(0,0))
# Cauchy regression using large values as starting values
reg5 <- cauchylm(c(y,y0),c(x,x0),c(10,10))



#Shows that EM algorithm converged to different modes
data.frame(
ols_start=sum(dcauchy(c(y,y0)-reg1$beta[1]-reg1$beta[2]*c(x,x0),0,reg1$sigma,log=TRUE)),
    true_start=sum(dcauchy(c(y,y0)-reg2$beta[1]-reg2$beta[2]*c(x,x0),0,reg2$sigma,log=TRUE)),
rnorm_start=sum(dcauchy(c(y,y0)-reg3$beta[1]-reg3$beta[2]*c(x,x0),0,reg3$sigma,log=TRUE)),
    zero_start=sum(dcauchy(c(y,y0)-reg4$beta[1]-reg4$beta[2]*c(x,x0),0,reg4$sigma,log=TRUE)),
large_start=sum(dcauchy(c(y,y0)-reg5$beta[1]-reg5$beta[2]*c(x,x0),0,reg5$sigma,log=TRUE))
)

#plot of the fitted lines
matplot(c(x,x0),cbind(
reg1$ols_reg$fitted,
reg2$ols_reg$fitted,
reg3$ols_reg$fitted,
reg4$ols_reg$fitted,
reg5$ols_reg$fitted),type='l',ylab="Y",xlab="X",main="Cauchy Regression")
points(c(x,x0),c(y,y0))
legend(x0-1,2,c("OLS","TRUE","RAND","ZERO","LARGE"),col=1:5,lty=1:5)
$\endgroup$
4
  • $\begingroup$ @probilityislogic: which mean the optimisation problem is not convex -so you may have multiple modes. One question (this is actually something I was trying to elucidate). Reading On the uni-modality of the likelihood for the Cauchy distribution (Copas, 1975). Biometrika 62, 701-704. I was under the impression that, except for a pathological case of bi-modality, the solution to this problem is unique. I was wondering whether you could comment on this (sorry if this sounds off base but this question is what got me interested in the Cauchy M in the first place). $\endgroup$
    – user603
    Apr 29 '13 at 15:53
  • $\begingroup$ This paper only deals with the case of no covariates. If you think in terms of IRLS the negative hessian has the form $\sum_iw_ix_ix_i^T$. For Cauchy the weights are positive for "good" data with small residuals and negative for "bad" data with large residuals. The weights are bounded though, so as long as there aren't too many outliers you should be fine. $\endgroup$ Apr 30 '13 at 8:16
  • $\begingroup$ so the negative hessian is PSD....and the optimization problem is convex? I don't understand your last sentence: did you mean 'fine' in wrt the convexity or some other aspect of the problem? $\endgroup$
    – user603
    Apr 30 '13 at 8:22
  • 1
    $\begingroup$ I'll post some R code giving an example. The negative hessian is PSD at the mode(s) but there may be more than one mode. Depending on where you start the algorithm you get different stationary points. $\endgroup$ May 1 '13 at 7:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.