Modelo cuando errores toman una distribución de Cauchy

Question

Es mi entendimiento de que la suma de los cuadrados de los errores (ESS) sirve como un estimador de máxima verosimilitud cuando un modelo de errores están distribuidos normalmente. (Es decir, si usted encontrar los parámetros del modelo que minimice la ESS, también maximizar la probabilidad.) Sin embargo, la distribución de error de mi modelo se parece mucho más como una distribución de Cauchy. Sería minimizar la ESS, el resultado de la máxima verosimilitud conjunto de parámetros para mi modelo? Si no, ¿qué estadística debo mirar?

Perdóname si esto no tiene ningún sentido, o me estoy perdiendo algo simple. Por favor, siéntase libre de enlazar a las fuentes que me puede ayudar a entender los conceptos básicos un poco mejor. Gracias!

Answer 1

Las estimaciones de mínimos cuadrados de los coeficientes de regresión son sólo igual a la máxima de estimaciones de probabilidad cuando los errores tienen una distribución normal (ver aquí para la prueba).

Si realmente quería máximo de estimaciones de probabilidad de los parámetros de regresión con Cauchy errores, basta con mirar a la probabilidad: $$L(\beta,\sigma)=\prod_{i=1}^n {\frac{1}{\pi\sigma\left(1+\left(\frac{y_i-\beta^\mathrm{T}x_i}{\sigma}\right)^2\right)}}$$

($y_i$$i$th observación, $x_i$ el vector de factores, $\sigma$ el parámetro de escala, & $\beta$ el vector de coeficientes.) No hay suficientes estadística de dimensión menor que el conjunto de datos completo, así que no es tan fácil de aprovechar al máximo, a pesar de que probablemente hay un método mejor que la fuerza bruta. Pero sin algunos teóricos de la motivación para la asunción de Cauchy errores, sólo se puede decir que tienen un poco de grasa de cola de la distribución. En esta situación de una forma u otra de regresión robusta podría ser vale la pena considerar.

Tenga en cuenta que el mínimo de plazas enfoque no es la cosa peor que usted podría utilizar incluso así. Siempre que la varianza es constante (y finito, que no lo es para la de Cauchy) todavía proporciona estimaciones coherentes, incluso el mejor lineales estimaciones imparciales, aunque tendrías que tomar intervalos de confianza con una pizca de sal.

Answer 2

El Jeffreys (posterior) de distribución es muy bueno para hacer inferencia en un modelo de regresión lineal de la ubicación de la escala de los errores. La inferencia basado en la Jeffreys distribución logra muy buena frecuentista propiedades: proporciona intervalos de confianza cuya cobertura es cercana a la nominal de la cobertura, incluso para tamaños de muestra pequeños.

Deje $\Phi$ ser cualquier derivable cdf y $\phi=\Phi'$ el correspondiente pdf. Considere el modelo de regresión lineal $Y=X\beta+\sigma\epsilon$$\epsilon_i \sim_{\text{iid}} \mathrm{d}\Phi$.

El Jeffreys distribución posterior es dada, a una constante de proporcionalidad, por: $$ \pi(\beta \sigma \mid y) \propto \frac{1}{\sigma^{n+1}} \prod_{i=1}^n \phi\left(\frac{y_i-\mu_i}{\sigma}\right) =: f(\beta \sigma \mid y) $$ donde $\mu_i=X_i \beta$ es el valor esperado de $y_i$.

El problema, dado un conjunto $A$ $q$- dimensional espacio de parámetros $\Theta$ ($q=p+1$desde los parámetros son $\beta_1$, $\ldots$, $\beta_p$ y $\sigma$), es evaluar $$ \int_A \pi(\beta \sigma \mid y) \mathrm{d}\beta \mathrm{d}\sigma = \frac{\int_A f(\beta \sigma \mid y)\mathrm{d}\beta \mathrm{d}\sigma}{\int_\Theta f(\beta \sigma \mid y)\mathrm{d}\beta\mathrm{d}\sigma} \aprox ? $$

Change of variables

It is possible to transform this integral to an integral on ${[0,1]}^q$ de la siguiente manera. El punto clave es el hecho de que $$ \frac{1}{\sigma^{q+1}}\prod_{i=1}^q\phi\left(\frac{y_i-\mu_i}{\sigma}\right) $$ es decir, hasta una constante de proporcionalidad no dependiendo de la $(\beta,\sigma)$, el Jacobiano de la función $$ F\colon(\beta\sigma)\mapsto \left(\Phi\left(\frac{y_1-\mu_1}{\sigma}\right), \ldots, \Phi\left(\frac{y_q-\mu_q}{\sigma}\right)\right) \in {[0,1]}^q. $$ No he tratado de demostrar este punto, pero se puede comprobar numéricamente, para una regresión lineal simple, por ejemplo:

x <- 1:10; y <- rcauchy(length(x), 2+x)
X <- model.matrix(~x)
Phi <- pcauchy; phi <- dcauchy # use any cdf and pdf you want
F <- function(betasigma){
  sapply(1:3, function(i) Phi((y[i]-X[i,]%*%betasigma[1:2])/betasigma[3]))
}
library(pracma) # provides the jacobian() function
f <- function(betasigma){
  prod(sapply(1:3, function(i){ phi((y[i]-X[i,]%*%betasigma[1:2])/betasigma[3]) })) / betasigma[3]^4
}
# look, the ratio is always the same:
det(jacobian(F, c(1,1,1)))/f(c(1,1,1))
## [1] -19.01263
det(jacobian(F, c(1,2,1)))/f(c(1,2,1))
## [1] -19.01263
det(jacobian(F, c(2,2,2)))/f(c(2,2,2))
## [1] -19.01263

Por lo tanto, $$ \begin{align} \int_A f(\beta, \sigma \mid y)d\beta d\sigma & \propto \int_A \bigl|\det J_F(\mu,\sigma)\bigr| \frac{1}{\sigma^{n-q}} \prod_{i=q+1}^n \phi\left(\frac{y_i-\mu_i}{\sigma}\right) \mathrm{d}\beta\mathrm{d}\sigma \\ & = \int_{F(A)} g\bigl(F^{-1}(u_1, \ldots, u_q)\bigr)\mathrm{d}u_1\ldots\mathrm{d}u_q \end{align} $$ donde $g(\beta,\sigma)=\frac{1}{\sigma^{q+1}} \prod_{i=1}^q \phi\left(\frac{y_i-\mu_i}{\sigma}\right)$.

No es difícil obtener la inversa de a $F$: $$ F^{-1}(u_1, \ldots, u_q) = {(\beta\sigma)}' = {(H, H)}^{-1}H'y_{1:q} $$ donde la matriz $H$$H=\left[ X_{1:q}, {\bigl(\Phi^{-1}(u_i)\bigr)}_{i\in(1:q)}\right]$.

Tenga en cuenta que $F^{-1}(u_1, \ldots, u_q)$ rendimientos $\sigma<0$ para algunos valores de la $u_i$. De hecho, si $F^{-1}\bigl(\Phi(z_1), \ldots, \Phi(z_q))={(\beta,\sigma)}'$,$F^{-1}\bigl(\Phi(-z_1), \ldots, \Phi(-z_q))={(\beta,-\sigma)}'$, por lo tanto, el conjunto de $u_i$'s para que $\sigma>0$ ha Lesbegue medida $1/2$.

De hecho, la Jeffreys distribución para una ubicación en la escala de regresión lineal es la misma que la fiduciaria de distribución. El método que se presenta es un caso particular del método general dada en el papel Computacional de los problemas de la generalizada fiduciales de inferencia por Hannig & al.. Pero hay un alto simplificación del método general en el caso de una ubicación a escala de regresión lineal (podemos tomar $K=1$ con las notaciones de la de papel, pero no voy a desarrollar este punto).

Algoritmo

El Jeffreys siguiente función devuelve una aproximación de la Jeffreys de distribución para el modelo de regresión lineal cuando los errores siguen una distribución para Estudiantes con grados de libertad df, a ser establecido por el usuario. Para df=Inf (por defecto), este es el Gaussiano de regresión lineal; por df=1 este es el Cauchy de regresión lineal. En el caso Gaussiano df=Inf, podemos comparar los resultados con la exacta Jeffreys distribución, que es conocido y primaria. Además, la inferencia basada en el Jeffreys distribución en el caso Gaussiano es el mismo que el habitual de mínimos cuadrados de inferencia (como veremos en los ejemplos).

De forma predeterminada, el X de la matriz es la matriz de la intersección modelo sólo y~1.

La aproximación se obtiene por una de Riemann-como la integración en ${[0,1]}^q$ mediante el uso de un uniforme de la partición en hypercubes. La partición está controlado por el argumento L, dando el número de centros de la hypercubes en cada coordenada (por lo tanto, hay $L^q$ hypercubes).

#' parameters: y (sample), X (model matrix), L (number of points per coordinate)
Jeffreys <- function(y, X=as.matrix(rep(1,length(y))), L=10, df=Inf){
  qdistr <- function(x, ...) qt(x, df=df, ...)
  ddistr <- function(x, ...) dt(x, df=df, ...)
  n <- nrow(X)
  q <- ncol(X)+1
  # centers of hypercubes (volume 1/L^p)
  centers <- as.matrix(do.call(expand.grid, rep(list(seq(0, 1, length.out=L+1)[-1] - 1/(2*L)), q)))
  # remove centers having equal coordinates (H'H is not invertible)
  centers <- centers[apply(centers, 1, function(row) length(unique(row))>1),]
  # outputs
  M <- (L^q-L)/2 # number of centers yielding sigma>0
  J <-  numeric(M)
  Theta <- array(0, c(M, q))
  # algorithm
  I <- 1:q
  yI <- y[I]; ymI <- y[-I]
  XI <- X[I,]; XmI <- X[-I,]
  counter <- 0
  for(m in 1:nrow(centers)){
    H <- unname(cbind(XI, qdistr(centers[m,])))
    theta <- solve(crossprod(H))%*%t(H)%*%yI
    if(theta[q]>0){ # sigma>0
      counter <- counter+1
      J[counter] <- sum(ddistr((ymI-XmI%*%head(theta,-1))/theta[q], log=TRUE)) - (n-q)*log(theta[q])  
      Theta[counter,] <- theta
    }
  }
  J <- exp(J)
  return(list(Beta=Theta[,-q], sigma=Theta[,q], W=J/sum(J)))
}

La función devuelve los valores de $(\beta,\sigma)$ correspondiente a cada hipercubo en el centro de la partición de ${[0,1]}^q$. También calcula los valores del integrando evaluados en cada centro en el vector J, y devuelve el vector normalizado de pesos W=J/sum(J). Vamos a ver cómo tratar con estas salidas en algunos ejemplos.

Primer ejemplo: Gaussiano de la muestra

Vamos a probarlo por un yo.yo.d. De gauss de la muestra $y_i \sim_{\text{i.i.d.}} {\cal N}(\mu, \sigma^2)$:

set.seed(666)
n <- 4
y <- rnorm(n)
results <- Jeffreys(y, L=100)
Mu <- results$Beta; Sigma <- results$sigma; W <- results$W

Ahora podemos tratar Mu y Sigma como si fueran ponderado de las muestras de la Jeffreys distribución, con pesos W.

El teórico media es la media de la muestra, y nuestra aproximación es bastante buena:

sum(W*Mu); mean(y)
## [1] 1.109794
## [1] 1.110175

Podemos obtener el aproximado Jeffreys cdf con la ewcdf función (ponderado empírica cdf) de la spatstat paquete, y comparar con la teórica. Nuestra aproximación es bastante perfecto:

### approximate Jeffreys distribution of µ ###
F_mu <- spatstat::ewcdf(Mu, weights=W)
curve(F_mu, from=0, to=2.5, xlab="mu", ylim=c(0,1), col="blue", lwd=2)
### exact Jeffreys distribution ###
mean_y <- mean(y); sd_y <- sd(y)
curve(pt((x-mean_y)/(sd_y/sqrt(n)), df=n-1), add=TRUE, col="red", lwd=4, lty="dashed")

Podemos obtener los intervalos de confianza mediante la aplicación de la quantile función de la ponderación de la cdf F_mu. Ellos son, en teoría, el mismo que los dolores de los intervalos de confianza en Gaussiano de regresión lineal, y realmente llegamos muy cerca de los resultados:

quantile(F_mu, c(2.5,97.5)/100)
##       2.5%      97.5% 
## -0.7143989  2.9309891
confint(lm(y~1)) # theoretically the same
##                  2.5 %  97.5 %
## (Intercept) -0.7121603 2.93251

El mismo para $\sigma$ (conocer el inverso de distribución Gamma de $\sigma^2$):

F_sigma <- spatstat::ewcdf(Sigma,W)
curve(F_sigma, from=0, to=2.5, xlab="sigma", ylim=c(0,1), col="blue", lwd=2)
curve(1-pgamma(1/x^2, (n - 1)/2, (n - 1) * sd_y^2/2), add=TRUE, col="red", lwd=4, lty="dashed")

Segundo ejemplo: Cauchy de la muestra

Ahora vamos a probar un yo.yo.d. Cauchy de la muestra con el tamaño de la muestra $n=200$.

set.seed(666)
n <- 200
y <- rcauchy(n)
results <- Jeffreys(y, L=100, df=1)
Mu <- results$Beta; Sigma <- results$sigma; W <- results$W

Desde $n=200$ no es un tamaño pequeño de la muestra, el Jeffreys significa que están cerca del máximo de estimaciones de probabilidad:

sum(W*Mu); sum(W*Sigma)
## [1] -0.01490355
## [1] 0.9081371
MASS::fitdistr(y, "cauchy")
##     location        scale   
##   -0.01345121    0.89958785 
##  ( 0.09185580) ( 0.08874509)

El MASS::rlm estimaciones no son tan de cerca:

rlmfit <- MASS::rlm(y~1)
rlmfit$coefficients
    ## (Intercept) 
    ##  -0.1160915
    rlmfit$s # rlm estimate of sigma
## [1] 1.338744

Jeffreys intervalos de confianza están cerca de la ML asintótica de los intervalos de confianza:

F_mu <- spatstat::ewcdf(Mu, weights=W); F_sigma <- spatstat::ewcdf(Sigma,W)
quantile(F_mu, c(2.5,97.5)/100)
##       2.5%      97.5% 
## -0.1971883  0.1707172
quantile(F_sigma, c(2.5,50,97.5)/100)
##      2.5%       50%     97.5% 
## 0.7471966 0.9055118 1.1027395
confint(MASS::fitdistr(y, "cauchy"))
##               2.5 %    97.5 %
## location -0.1934853 0.1665829
## scale     0.7256507 1.0735250

Tercer ejemplo : Gaussiano de regresión lineal simple

Niza:

f <- function(x) 4+2*x
set.seed(666)
n <- 20
x <- seq_len(n) # covariates
y <- f(x)+rnorm(n)
# run algorithm
results <- Jeffreys(y, X=model.matrix(~x), L=60)
# outputs
W <- results$W; Beta0 <- results$Beta[,1]; Beta1 <- results$Beta[,2]
sum(W*Beta0); sum(W*Beta1) # Jeffreys means
## [1] 4.172721
## [1] 1.983503
coef(lm(y~x)) # theoretically the same
## (Intercept)           x 
##    4.179859    1.983008
F_Beta0 <- spatstat::ewcdf(Beta0, weights=W); F_Beta1 <- spatstat::ewcdf(Beta1, weights=W)
quantile(F_Beta0, c(2.5,97.5)/100); quantile(F_Beta1, c(2.5,97.5)/100)
##     2.5%    97.5% 
## 2.883869 5.499620
##     2.5%    97.5% 
## 1.872328 2.095764
confint(lm(y~x)) # theoretically the same
##                2.5 %   97.5 %
## (Intercept) 2.857903 5.501815
## x           1.872653 2.093362

Cuarto ejemplo : Cauchy de regresión lineal simple

set.seed(666)
y <- f(x)+rcauchy(n)
# run algorithm
results <- Jeffreys(y, X=model.matrix(~x), L=60, df=1)
# outputs
W <- results$W; Beta0 <- results$Beta[,1]; Beta1 <- results$Beta[,2]; Sigma <- results$sigma
# Jeffreys means
sum(W*Beta0); sum(W*Beta1); sum(W*Sigma)
## [1] 4.157664
## [1] 1.997121
## [1] 0.685825

Mientras que $n=20$ no es grande, el ML estimaciones de los parámetros de regresión son cerca de sus Jeffreys medios, pero no tan cerca para $\sigma$:

X <- model.matrix(~x)
likelihood <- function(y, beta0, beta1, sigma){
  prod(dcauchy((y-X%*%c(beta0,beta1))/sigma)/sigma)
}
(ML <- MASS::fitdistr(y, likelihood, list(beta0=sum(W*Beta0), beta1=sum(W*Beta1), sigma=1)))
##      beta0        beta1        sigma   
##   4.20188590   1.99239112   0.60087397 
##  (0.54295228) (0.04433536) (0.18660186)

El Jeffreys intervalos de confianza están cerca de la ML intervalos de confianza, a excepción de $\sigma$:

confint(ML)
##          2.5 %    97.5 %
## beta0 3.137719 5.2660528
## beta1 1.905495 2.0792868
## sigma 0.235141 0.9666069
F_Beta0 <- spatstat::ewcdf(Beta0, weights=W); F_Beta1 <- spatstat::ewcdf(Beta1, weights=W)
quantile(F_Beta0, c(2.5,97.5)/100); quantile(F_Beta1, c(2.5,97.5)/100)
##     2.5%    97.5% 
## 3.098442 5.167351
##     2.5%    97.5% 
## 1.913328 2.089146
F_sigma <- spatstat::ewcdf(Sigma,W); quantile(F_sigma, c(2.5,50,97.5)/100)
##      2.5%       50%     97.5% 
## 0.3418978 0.6491162 1.2464163

El MASS::rlm las estimaciones de los parámetros de regresión son más bien cerca de sus Jeffreys significa también:

rlmfit <- MASS::rlm(y~x)
rlmfit$coefficients
    ## (Intercept)           x 
    ##    3.945603    2.042590
    rlmfit$s # rlm estimate of sigma
## [1] 1.019974

Answer 3

GraphPad Prism puede hacer la regresión no lineal, asumiendo una distribución de Cauchy. Ese es nuestro método robusto. La matemática detalles están explicados en detalle en las páginas 11-14 de BMC Bioinformatics 2006, 7:123 doi:10.1186/1471-2105-7-123, la Detección de valores atípicos cuando los datos de adaptación con la regresión no lineal – un nuevo método basado en la solidez de regresión lineal y el falso tasa de descubrimiento

Answer 4

El uso de Cauchy errores NO ES un método robusto. Esto conduce a un modelo que puede capturar los valores extremos, pero si no hay valores atípicos, el resultado del modelo es muy restrictiva, ya que se supone que la distribución de los errores es pesado de cola con una cola de comportamiento. Debido a que la distribución de Cauchy es un caso especial de la distribución t para $\nu=1$, esto hace que una muy fuerte declaración acerca de cómo los errores se distribuyen. Una actuación más eficaz consiste en la utilización de un Estudiante $t$- $t(0,\sigma,\nu)$ donde $\nu$ son los grados de libertad y $\sigma$ es el parámetro de escala, que son desconocidos y deben ser estimada a partir de los datos.

El modelo es

$$y_j = h(x_j^{\top}\beta) + e_j,$$

donde $e_j\stackrel{ind}{\sim} t(0,\sigma,\nu)$, $j=1,\dots,n$, y $h$ es una función real.

La probabilidad de $(\beta,\sigma,\nu)$ es el dado por

$${\mathcal L}(\beta,\sigma,\nu) \propto \prod_{j=1}^n f(h(y_j- x_j^{\top}\beta);0,\sigma,\nu),$$

donde

$$f(z;\mu,\sigma,\nu) = \dfrac{1}{\sigma}\dfrac{\Gamma\left(\dfrac{\nu+1}{2}\right)}{\sqrt{\pi\nu}\Gamma\left(\dfrac{\nu}{2}\right)} \left[1+\dfrac{1}{\nu}\left(\dfrac{z-\mu}{\sigma}\right)^2\right]^{-\frac{\nu+1}{2}}.$$

Los estimadores de máxima verosimilitud puede ser obtenido utilizando métodos numéricos. Tenga en cuenta que esta estructura cubre el ajuste de datos, lineales y no lineales de regresión.

De La Wikipedia:

Estadísticas robustas son las estadísticas con un buen rendimiento para los datos extraídos de una amplia gama de distribuciones de probabilidad, especialmente para las distribuciones que no están normalmente distribuidos. Estadístico robusto se han desarrollado métodos para muchos problemas comunes, tales como la estimación de ubicación, escala y los parámetros de regresión. La motivación es la producción de métodos estadísticos que no son excesivamente afectada por los valores extremos. Otra motivación es proporcionar métodos con un buen rendimiento cuando hay pequeñas desviaciones de la paramétrica de las distribuciones. Por ejemplo, robusto métodos funcionan bien para las mezclas de dos distribuciones normales con diferentes estándar de las desviaciones, por ejemplo, de uno y tres; bajo este modelo, no de métodos robustos como una prueba de t de funcionar mal.