Ajuste de distribuciones personalizadas mediante MLE

Question

Mi pregunta está relacionada con el ajuste de distribuciones personalizadas en R pero creo que tiene suficiente elemento probabilístico como para quedarse en CV.

Tengo un interesante conjunto de datos que presenta las siguientes características:

• Masa grande a cero
• Masa considerable por debajo de un umbral que se ajusta muy bien a una distribución paramétrica sesgada a la derecha.
• Pequeña cantidad de masa en valores extremos
• Una serie de covariables que deben impulsar la variable de interés

Esperaba poder modelizar esto utilizando un enfoque de distribución cero-inflada, ampliamente explorado en la literatura. Esencialmente, la densidad es:

$$f_{Y}(y)=\begin{cases} \pi \quad\quad\quad\quad\,\,\,\,\,\,,\,\,y=0 \\ (1-\pi)f_X(y),\,\,y>0 \end{cases}$$

Esto es bastante fácil de instalar tal cual. Sin embargo, me gustaría que el parámetro de mezcla $\pi$ dependa de las covariables $Z$ mediante una regresión logística:

$$\text{logit}(\mathbb{E}[\pi\,|\,Z])=\beta Z$$ donde $\beta$ es un vector de coeficientes para las covariables. Además, debido a la naturaleza de cola extrema de mis datos, mi distribución $f_{X}(y)$ encaja mejor con un enfoque de valor extremo:

$$f_{X}(y)=\begin{cases} f_{A}(y;a,b) \quad\,\,\,\,\,\,\,,\,y\leq \mu \\ (1-F_{A}(\mu))\cdot\text{GPD}\bigg(y;\mu,\sigma=\frac{(1-F_{A}(\mu))}{f_{A}(\mu)},\xi\bigg),\,y>\mu \end{cases}$$ donde $\text{GPD}(y;\mu,\sigma,\xi)$ se refiere a la distribución de Pareto generalizada, que modela el exceso por encima de un determinado umbral $\mu$ y $f_{A}(y;a,b)$ es una distribución asimétrica dada con parámetros de escala y forma $a$ y $b$ respectivamente. La caracterización anterior garantiza que las densidades son continuas en $y=\mu$ (aunque no diferenciable) y que $f_{X}(y)$ se integra en 1.

Además, lo ideal sería que los parámetros de las distribuciones anteriores también dependieran de las covariables:

$$f_{A}(y;a,b,\beta Z)$$ $$\text{GPD}\bigg(y;\mu,\sigma=\frac{(1-F_{A}(\mu))}{f_{A}(\mu)},\xi,\beta Z\bigg)$$

Me doy cuenta de que la configuración anterior es bastante compleja, pero me preguntaba si hay una manera de derivar las estimaciones MLE de cada uno de los parámetros deseados mediante la maximización de la función de verosimilitud es decir, para obtener:

$$\hat{\xi}, \hat{a}, \hat{b}, \hat{\beta}$$

¿Hay alguna forma factible/ideal de hacerlo en R? Tanto en lo que respecta a mi problema específico como al ajuste de distribuciones personalizadas en general.

Answer 1

Esta respuesta presupone $\mu$ es conocido.

Una forma muy flexible de obtener MLE en R es utilizar STAN mediante rstan . STAN tiene fama de ser una herramienta MCMC, pero también puede estimar parámetros mediante inferencia variacional o MAP. Y usted es libre de no especificar los antecedentes .

En este caso, lo que estás haciendo es muy similar a su ejemplo de modelo de obstáculos . Aquí está el código STAN para ese ejemplo.

data {
  int<lower=0> N;
  int<lower=0> y[N];
}
parameters {
  real<lower=0, upper=1> theta;
  real<lower=0> lambda;
}
model {
  for (n in 1:N) {
    if (y[n] == 0)
      target += bernoulli_lpmf(1 | theta);
    else
      target += bernoulli_lpmf(0 | theta)
                  + poisson_lpmf(y[n] | lambda);
  }
}

Para adaptar esto a tu propio uso, podrías:

• Sustituir poisson_lpmf con la densidad logarítmica para su $f_A$ .
• Añade un tercer caso al if-else para que compruebe si se supera $\mu$ Como la carne de ese tercer caso, utilice el log pmf para su distribución de valor extremo de elección.
• Sustituir bernoulli_lpmf con categorical_lpmf y convertir el parámetro de probabilidad de la mezcla en un vector.
• Para incorporar covariables, puede añadir parámetros de regresión y hacer que todos los demás parámetros sean funciones de ellos. Puede ser útil utilizar categorical_logit_lpmf en lugar de categorical_lpmf .
• Truncar un componente de la mezcla a $\mu$ desde arriba y el otro en $\mu$ de abajo, en función de tu perspectiva sobre el dilema planteado por Jarle Tufto en los comentarios. Parece que se podrían obtener estimaciones MUY diferentes dependiendo de cómo se decida exactamente manejar esto. Una buena comprobación: generar un conjunto de datos falso a partir de los parámetros ajustados y asegurarse de que tiene la cantidad correcta en 0, la cantidad por encima de $\mu$ etc.

Una vez que tenga un archivo con el código STAN correcto, puede utilizar STAN con muchas cadenas de herramientas diferentes. Para utilizarlo con R, echa un vistazo a estos ejemplos . He simplificado una para obtener una MLE, utilizando rstan::optimizing en lugar de sampling :

install.packages("rstan")
library("rstan")
model = stan_model("Example1.stan")
fit = optimizing(model)

También hay algunos trucos para una optimización más rápida/mejor que podría ayudar en la práctica.

Answer 2

Mi STANswer es tan compleja que está pidiendo a gritos que algo salga mal. Aquí hay una manera más simple: hacer toda su inferencia condicional en los hechos (conocidos) de si cada dato es superior a 0 y si cada dato es superior a $\mu$ . En otras palabras, reducir los datos a:

• El conjunto de observaciones $S_1 \equiv \{y_i: 0<x_i<\mu\}$ .
• El conjunto de observaciones $S_2 \equiv \{y_i: \mu<x_i\}$ .
• El número de ceros $N_0$ .
• Sea $N_1 \equiv |S_1|$ y $N_2 \equiv |S_2|$ .

A continuación, maximizar, por separado:

• $(1-F_A(\mu))^{N_1} \prod_{y \in S_1} f_A(y) $ en relación con los parámetros de $f_A$ .
• $\prod_{y \in S_2} GPD(y)$ en relación con los parámetros de su promesticoducto bruto (GPD).
• $\pi^{N_0}(1-\pi)^{N_1 + N_2}$ por ejemplo $\pi$ .

Esto no parece hacer exactamente lo que pides porque:

• $\mu$ debe ser especificado por el usuario.
• el parámetro de escala GPD no se fija en $\frac{f_A(\mu)}{1-F_A(\mu)}$ . Esperemos que esa parte no sea esencial para la interpretabilidad. Si lo es, tal vez bastaría con arreglar $\sigma$ basándose en los resultados del punto 1, y luego optimizar el resto de parámetros. Entonces ya no es un MLE conjunto. No hay manera de desacoplar la optimización si usted está realmente decidido a eso.