Distribución normal con muchos valores cero

Question

Tengo aquí una función de densidad de probabilidad que representa la cantidad de horas dedicadas al estudio para un examen final de clase. No sorprendentemente, muchos de los estudiantes no estudian para el final.

Sin embargo, debido a que muchos estudiantes acumulan 0 horas estudiadas, parece que la curva se extiende hacia un dominio negativo, lo cual no es posible ya que no puedes estudiar por horas negativas.

Aquí algunas preguntas:

• ¿Cómo se puede interpretar los resultados para una curva como esta?
• ¿Se pueden calcular de manera confiable valores de probabilidad, puntaje z, desviación estándar, función de distribución acumulada, etc?
• ¿Existe alguna forma de transformar esto para que todos los valores del dominio sean mayores a 0?

Answer 1

Creo que sería mejor tratar esto como una mezcla de dos distribuciones en lugar de intentar aplicar las herramientas estándar de la teoría normal. Por lo tanto, voy a esbozar un poco sobre la distribución Gamma inflada en ceros, incluyendo el cálculo de sus dos primeros momentos, para darte una idea de cómo va esto. Podrías cambiar fácilmente la Gamma por una distribución continua diferente si así lo prefieres (por ejemplo, una distribución Beta escalada para estar en $[0,100]$). Estoy dispuesto a agregar actualizaciones más adelante si esto no es útil.

Sea $Z_1,\dots,Z_n\stackrel{\text{iid}}\sim\text{Bern}(\theta)$ y considera $X_1,\dots,X_n$ donde $$ X_i \vert Z_i \stackrel{\text{iid}}\sim \begin{cases}\Gamma(\alpha, \beta) & Z_i = 1 \\ 0 & Z_i = 0\end{cases} $$

Entonces cada $X_i$ es una mezcla de una masa puntual en $0$ con probabilidad $1-\theta$ y una $\Gamma(\alpha,\beta)$ con probabilidad $\theta$. Interpretamos esto como que $Z_i$ es una variable latente oculta que determina si el estudiante estudia o no, y luego $X_i$ es el valor observado.

Esto es un poco formal, pero lo mencionaré por completitud. $X_i$ no tiene una función de densidad de probabilidad en el sentido habitual porque no es ni discreta ni continua, pero si consideramos la medida $\nu = \lambda + \delta_0$, es decir, la medida de Lebesgue más una masa puntual en $0$, entonces $\nu(A) = 0 \implies P_X(A) = 0$ para cualquier conjunto medible $A$, por lo que podemos obtener una función de densidad de probabilidad $f_X := \frac{\text dP_X}{\text d\nu}$ con respecto a $\nu$.

Pero, ¿cómo luce esta función de densidad de probabilidad? Podemos calcular la función de distribución acumulada (CDF) $F$ usando algunas reglas de probabilidad condicional. $$ F(x) = P(X\leq x \cap Z = 0) + P(X\leq x \cap Z = 1) \\ = P(X\leq x \vert Z = 0)P(Z=0) + P(X\leq x \vert Z=1)P(Z=1) \\ = 1 \cdot (1 - \theta) + F_\Gamma(x; \alpha, \beta) \theta \\ = 1 - \theta + \theta F_\Gamma(x; \alpha, \beta) $$ donde $F_\Gamma$ denota la CDF de una distribución Gamma real.

Entonces queremos una función $f_X$ tal que $$ F(x) = \int_{[0, x]} f_X\,\text d\nu. $$ Nota que $$ \int_{[0, x]} f_X\,\text d\nu = \int_{\{0\}} f_X\,\text d\delta_0 + \int_{(0, x)} f_X\,\text d\lambda \\ = f_X(0) + \int_{(0, x)} f_X\,\text d\lambda $$ así que puedo tomar $$ f_X(x) = 1 - \theta + \theta f_\Gamma(x; \alpha, \beta). $$

Veamos que esta es una función de densidad de probabilidad válida: $$ \int_{[0,\infty)} f_X\,\text d\nu = 1 - \theta + \theta \int_0^\infty f_\Gamma \,\text d\lambda = 1 $$ por lo que esta es de hecho una función de densidad de probabilidad válida (con respecto a $\nu$).

---

Ahora calcularé los dos primeros momentos de $X_i$.

$$ E(X_i) = \int_{[0,\infty)} x f_X(x)\,\text d\nu(x) \\ = 0 \cdot \int_{\{0\}} f_X\,\text d\delta_0 + \int_{(0,\infty)} x f_X(x)\,\text d\lambda(x) \\ = 0 + \theta \int_0^\infty x f_\Gamma(x) \,\text d\lambda(x) = \frac{\theta\alpha}\beta := \mu < \infty. $$ Luego $$ E(X_i^2) = \int_{[0,\infty)} x^2 f_X(x)\,\text d\nu(x) \\ = 0 + \theta \int_0^\infty x^2 f_\Gamma(x)\,\text d\lambda(x) \\ = \frac{\theta\alpha(1 + \alpha)}{\beta^2}. $$ Esto significa $$ \sigma^2 := E(X_i^2) - \mu^2 < \infty. $$

Por fin he confirmado los siguientes hechos: tenemos una colección de variables aleatorias independientes $X_1, X_2,\dots$ con medias y varianzas finitas, por lo que podemos aplicar felizmente el TCL estándar para concluir

$$ \frac{\bar X_n - \mu_n}{\sqrt n} \stackrel{\text d}\to \mathcal N(0, \sigma^2). $$

Ahora, en cuanto a qué tan bueno es esto, probablemente querrás hacer algunas simulaciones. Además, no estoy diciendo que este sea realmente un buen modelo.

---

Comprobaré mi matemática con la siguiente simulación:

theta <- .76
a <- 5.4
b <- 1.2
n <- 1e6

set.seed(42)
z <- rbinom(n, 1, theta)
x <- numeric(n)
x[z==1] <- rgamma(sum(z), shape=a, rate=b)

hist(x, main="Simulaciones de una Gamma inflada en cero")

mean(x)
theta * a / b  # coincide

mean(x^2)
theta * a * (1 + a) / b^2  # coincide

También ten en cuenta que no estoy utilizando una KDE para mostrar la distribución como (parece que) estás haciendo tú. Por lo general, esas no son apropiadas para distribuciones que tienen una masa puntual como esta. Además, si estás usando una que pone una mini gaussiana en cada punto de datos, entonces se asume implícitamente que el soporte es todo $\mathbb R$, por lo que también puedes obtener probabilidad positiva en áreas imposibles como lo hiciste.

Si decides usar este modelo y quieres estimar los parámetros, el algoritmo de esperanza-maximización es el camino habitual a seguir.

En este caso, sin embargo, no hay duda acerca de a qué clase pertenece un $X_i$ en particular porque si $X_i = 0$ entonces $Z_i = 0$ casi seguramente. Por lo tanto, puedes hacer:

mean(x > 0)  # comparar con theta

mu.x <- mean(x[x > 0])
s2.x <- var(x[x > 0])
(b.hat <- mu.x / s2.x)
(a.hat <- mu.x^2 / s2.x)

y estos coinciden. Pero tengo un tamaño de muestra enorme y $\theta$ no está particularmente cerca de $0$ o $1$ aquí, por lo que no es impresionante ser tan preciso con este enfoque condicional.