Ajuste de la distribución t en R: parámetro de escala

¿Cómo ajusto los parámetros de una distribución t, es decir, los parámetros correspondientes a la 'media' y 'desviación estándar' de una distribución normal? Supongo que se llaman 'media' y 'escala / grados de libertad' para una distribución t?

El siguiente código a menudo produce errores de "error de optimización".

library(MASS)
fitdistr(x, "t")

¿Tengo que escalar x primero o convertirlo en probabilidades? ¿Cómo mejor hacer eso?

fitdistrutiliza técnicas de máxima verosimilitud y optimización para encontrar parámetros de una distribución dada. A veces, especialmente para la distribución t, como notó @ user12719, la optimización en la forma:

fitdistr(x, "t")

falla con un error.

En este caso, debe ayudar al optimizador proporcionando el punto de partida y el límite inferior para comenzar a buscar parámetros óptimos:

fitdistr(x, "t", start = list(m=mean(x),s=sd(x), df=3), lower=c(-1, 0.001,1))

Tenga en cuenta que df=3es su mejor estimación de lo que dfpodría ser un "óptimo" . Después de proporcionar esta información adicional, su error desaparecerá.

Un par de extractos para ayudarlo a comprender mejor la mecánica interna de fitdistr:

Para las distribuciones Normal, log-Normal, geométrica, exponencial y de Poisson, se usan los MLE de forma cerrada (y los errores estándar exactos), y startno se deben suministrar.

...

Para las siguientes distribuciones con nombre, se calcularán valores iniciales razonables si startse omiten o solo se especifican parcialmente: "cauchy", "gamma", "logístico", "binomial negativo" (parametrizado por mu y tamaño), "t" y "weibull ". Tenga en cuenta que estos valores iniciales pueden no ser lo suficientemente buenos si el ajuste es deficiente: en particular, no son resistentes a los valores atípicos a menos que la distribución ajustada sea de cola larga.

MASS, el libro (4a edición, página 110) desaconseja tratar de estimar , el parámetro de grados de libertad en la distribución con máxima probabilidad (con algunas referencias bibliográficas: Lange et al. (1989), "Modelo estadístico robusto Usando la distribución t", JASA , 84 , 408 , y Fernandez & Steel (1999), "multivariante Student t modelos de regresión: trampas y la inferencia", Biometrika , 86 , 1 ).$\nu$$t$

La razón es que la función de probabilidad para basada en la función de densidad t, puede ser ilimitada y, en esos casos, no dará un máximo bien definido. Veamos un ejemplo artificial donde se conoce la ubicación y la escala (como la distribución estándar ) y solo se desconocen los grados de libertad. A continuación se muestra un código R, que simula algunos datos, define la función de probabilidad de registro y la traza:$\nu$$t$

set.seed(1234)
n <- 10
x <- rt(n,  df=2.5)

make_loglik  <-  function(x)
    Vectorize( function(nu) sum(dt(x, df=nu,  log=TRUE)) )

loglik  <-  make_loglik(x)
plot(loglik,  from=1,  to=100,  main="loglikelihood function for df     parameter", xlab="degrees of freedom")
abline(v=2.5,  col="red2")

Si juega con este código, puede encontrar algunos casos en los que hay un máximo bien definido, especialmente cuando el tamaño de la muestra es grande. ¿Pero el estimador de máxima verosimilitud es bueno?$n$

Probemos algunas simulaciones:

t_nu_mle  <-  function(x) {
    loglik  <-  make_loglik(x)
    res  <-  optimize(loglik, interval=c(0.01, 200), maximum=TRUE)$maximum
    res   
}

nus  <-  replicate(1000, {x <- rt(10, df=2.5)
    t_nu_mle(x) }, simplify=TRUE)

> mean(nus)
[1] 45.20767
> sd(nus)
[1] 78.77813

Mostrar la estimación es muy inestable (mirando el histograma, una parte considerable de los valores estimados se encuentra en el límite superior dado para optimizar 200).

Repetir con un tamaño de muestra más grande:

nus  <-  replicate(1000, {x <- rt(50, df=2.5)
    t_nu_mle(x) }, simplify=TRUE)
> mean(nus)
[1] 4.342724
> sd(nus)
[1] 14.40137

lo cual es mucho mejor, pero la media todavía está muy por encima del verdadero valor de 2.5.

Entonces recuerde que esta es una versión simplificada del problema real donde los parámetros de ubicación y escala también deben estimarse.

Si la razón de usar la distribución es "solidificar", entonces estimar partir de los datos puede destruir la robustez.$t$$\nu$

En la ayuda para fitdistr está este ejemplo:

fitdistr(x2, "t", df = 9)

indicando que solo necesita un valor para df. Pero eso supone estandarización.

Para un mayor control, también muestran

mydt <- function(x, m, s, df) dt((x-m)/s, df)/s
fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, lower = c(-Inf, 0))

donde los parámetros serían m = media, s = desviación estándar, df = grados de libertad