¿Puedo reconstruir una distribución normal a partir del tamaño de la muestra, y los valores mínimo y máximo? Puedo usar el punto medio como valor aproximado de la media

Question

Sé que esto puede ser un poco irregular, estadísticamente hablando, pero este es mi problema.

Tengo muchos datos de rango, es decir, el mínimo, máximo y tamaño de la muestra de una variable. Para algunos de estos datos también tengo una media, pero no muchos. Quiero comparar estos rangos entre sí para cuantificar la variabilidad de cada rango, y también comparar las medias. Tengo una buena razón para suponer que la distribución es simétrica en torno a la media, y que los datos tendrán una distribución gaussiana. Por esta razón estoy pensando en justificar el uso del punto medio de la distribución como un proxy de la media, cuando esta esté ausente.

Lo que quiero hacer es reconstruir una distribución para cada rango, y luego usar eso para proporcionar una desviación estándar o error estándar para esa distribución. La única información que tengo es el máximo y el mínimo observado de una muestra, y el punto medio como un proxy de la media.

De esta manera espero poder calcular medias ponderadas para cada grupo, y también determinar el coeficiente de variación para cada grupo también, basado en los datos de rango que tengo y mis suposiciones (de una distribución simétrica y normal).

Planeo utilizar R para hacer esto, así que cualquier ayuda con código sería apreciada también.

Answer 1

La función de distribución acumulativa conjunta para el mínimo $x_{(1)}$ y máximo $x_{(n)}$ de una muestra de $n$ de una distribución gaussiana con media $\mu$ y desviación estándar $\sigma$ es

$$ F(x_{(1)},x_{(n)};\mu,\sigma) = \Pr(X_{(1)}x_{(1)}, X_{(n)}

donde $\Phi(\cdot)$ es la función de distribución acumulativa gaussiana estándar. La diferenciación con respecto a $x_{(1)}$ y $x_{(n)}$ da como resultado la función de densidad de probabilidad conjunta

$$f(x_{(1)},x_{(n)};\mu,\sigma) =\\ n(n-1)\left[\Phi\left(\tfrac{x_{(n)}-\mu}{\sigma}\right) - \Phi\left(\tfrac{x_{(1)}-\mu}{\sigma}\right)\right]^{n-2}\cdot\phi\left(\tfrac{x_{(n)}-\mu}{\sigma}\right)\cdot\phi\left(\tfrac{x_{(1)}-\mu}{\sigma}\right)\cdot\tfrac{1}{\sigma^2}$$

donde $\phi(\cdot)$ es la función de densidad de probabilidad gaussiana estándar. Tomando el logaritmo y eliminando términos que no contienen parámetros se obtiene la función de log-verosimilitud

$$\ell(\mu,\sigma;x_{(1)},x_{(n)}) =\\ (n-2)\log\left[\Phi\left(\tfrac{x_{(n)}-\mu}{\sigma}\right) - \Phi\left(\tfrac{x_{(1)}-\mu}{\sigma}\right)\right] + \log\phi\left(\tfrac{x_{(n)}-\mu}{\sigma}\right) + \log\phi\left(\tfrac{x_{(1)}-\mu}{\sigma}\right) - 2\log\sigma$$

Esto no parece muy manejable, pero es fácil ver que se maximiza independientemente del valor de $\sigma$ al establecer $\mu=\hat\mu=\frac{x_{(n)}+x_{(1)}}{2}$, es decir, el punto medio, la primera término se maximiza cuando el argumento de una función de distribución acumulativa es el negativo del argumento de la otra; los segundo y tercer términos representan la verosimilitud conjunta de dos variables normales independientes.

Sustituir $\hat\mu$ en la log-verosimilitud y escribir $r=x_{(n)}-x_{(1)}$ da $$\ell(\sigma;x_{(1)},x_{(n)},\hat\mu)=(n-2)\log\left[1 - 2\Phi\left(\tfrac{-r}{2\sigma}\right)\right] - \frac{r^2}{4\sigma^2} -2\log{\sigma}$$

Esta expresión debe maximizarse numéricamente (por ejemplo, con optimize de paquete stat de R) para encontrar $\hat\sigma$. (Resulta que $\hat\sigma=k(n)\cdot r$, donde $k$ es una constante que depende solo de $n$—quizás alguien más matemáticamente hábil que yo podría mostrar por qué.)

Los estimados no sirven sin una medida de precisión acompañante. La información de Fisher observada se puede evaluar numéricamente (por ejemplo, con hessian del paquete numDeriv de R) y usar para calcular errores estándar aproximados:

$$I(\mu)=-\left.\frac{\partial^2{\ell(\mu;\hat\sigma)}}{(\partial\mu)^2}\right|_{\mu=\hat\mu}$$ $$I(\sigma)=-\left.\frac{\partial^2{\ell(\sigma;\hat\mu)}}{(\partial\sigma)^2}\right|_{\sigma=\hat\sigma}$$

Sería interesante comparar los estimados de verosimilitud y de método de momentos para $\sigma$ en términos de sesgo (¿es el MLE consistente?), varianza y error cuadrático medio. También está el problema de la estimación para aquellos grupos donde la media de la muestra se conoce además del mínimo y el máximo.

Answer 2

Necesitas relacionar el rango con la desviación estándar / varianza. Deja que $\mu$ sea la media, $\sigma$ la desviación estándar y $R=x_{(n)} - x_{(1)}$ sea el rango. Entonces, para la distribución normal tenemos que el $99.7$% de la masa de probabilidad se encuentra dentro de 3 desviaciones estándar de la media. Esto, como regla práctica, significa que con una probabilidad muy alta,

$$\mu + 3\sigma \approx x_{(n)}$$ y

$$\mu - 3\sigma \approx x_{(1)}$$

Restando el segundo del primero obtenemos

$$6\sigma \approx x_{(n)} - x_{(1)}= R$$ (esto, por cierto, es de donde viene la metodología de aseguramiento de calidad "seis sigma" en la industria). Luego puedes obtener una estimación para la desviación estándar por $$\hat \sigma = \frac 16 \Big(\bar x_{(n)} - \bar x_{(1)}\Big)$$ donde la barra denota promedios. Esto es cuando asumes que todas las submuestras provienen de la misma distribución (mencionaste tener rangos esperados). Si cada muestra es una normal diferente, con media y varianza diferentes, entonces puedes utilizar la fórmula para cada muestra, pero la incertidumbre / posible inexactitud en el valor estimado de la desviación estándar será mucho mayor.

Tener un valor para la media y para la desviación estándar caracteriza completamente la distribución normal.

Answer 3

Es fácil obtener la función de distribución de la máxima de la distribución normal (ver "P.max.norm" en el código). A partir de ella (con algo de cálculo) puedes obtener la función cuantil (ver "Q.max.norm").

Utilizando "Q.max.norm" y "Q.min.norm" puedes obtener la mediana del rango que está relacionada con N. Usando la idea presentada por Alecos Papadopoulos (en la respuesta anterior) puedes calcular la desviación estándar.

Intenta esto:

N = 100000    # el tamaño de la muestra

# Función de probabilidad dada q y N
P.max.norm <- function(q, N=1, mean=0, sd=1){
    pnorm(q,mean,sd)^N
} 
# Funciones cuantil dadas p y N
Q.max.norm <- function(p, N=1, mean=0, sd=1){
    qnorm(p^(1/N),mean,sd)
} 
Q.min.norm <- function(p, N=1, mean=0, sd=1){
    mean-(Q.max.norm(p, N=N, mean=mean, sd=sd)-mean)
} 

### vamos a probarlo (toma algo de tiempo)
Q.max.norm(0.5, N=N)  # La mediana en la máxima
Q.min.norm(0.5, N=N)  # La mediana en la mínima

iter = 100
median(replicate(iter, max(rnorm(N))))
median(replicate(iter, min(rnorm(N))))
# está bastante bien

### Vamos a intentar obtener estimaciones
true_mean = -3
true_sd = 2
N = 100000

x = rnorm(N, true_mean, true_sd)  # simulación
x.vec = range(x)                  # observaciones

# estimación
est_mean = mean(x.vec)
est_sd = diff(x.vec)/(Q.max.norm(0.5, N=N)-Q.min.norm(0.5, N=N))

c(true_mean, true_sd)
c(est_mean, est_sd)

# Bastante bueno, pero solo para N grande
# -3  2
# -3.252606  1.981593