¿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 para representar la media

Sé que esto podría ser un poco complicado, estadísticamente, pero este es mi problema.

Tengo muchos datos de rango, es decir, el tamaño mínimo, máximo y de 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 para comparar las medias. Tengo una buena razón para suponer que la distribución es simétrica alrededor de la media, y que los datos tendrán una distribución gaussiana. Por esta razón, creo que puedo justificar el uso del punto medio de la distribución como proxy de la media, cuando está ausente.

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

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

Planeo usar R para hacer esto, por lo que cualquier ayuda de código también sería apreciada.

r normal-distribution estimation missing-data order-statistics green_thinlake
fuente

Me preguntaba por qué dice que tiene datos para valores mínimos, máximos y máximos; luego, que tiene información sobre el mínimo y máximo esperado. ¿Cuál es - observado o esperado?

Scortchi - Restablece a Monica

Lo siento, ese es mi error. Se observan los datos máximos y mínimos (medidos a partir de objetos de la vida real). He modificado la publicación.

green_thinlake

Respuestas:

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

F (x_{(1)}, x_{(n)}; μ, σ) = Pr (X_{(1)} < x_{(1)}, X_{(n)} < x_{(n)}) = Pr (X_{(n)} < x_{(n)}) - Pr (X_{(1)} > x_{(1)}, X_{(n)} < x_{(n)} = Φ {(\frac{x_{(n)} - μ}{σ})}^{n} - {[Φ (\frac{x_{(n)} - μ}{σ}) - Φ (\frac{x_{(1)} - μ}{σ})]}^{n}

$F(x_{(1)},x_{(n)};\mu,\sigma) = \Pr(X_{(1)}<x_{(1)}, X_{(n)}<x_{(n)})\\ =\Pr( X_{(n)}<x_{(n)}) - \Pr(X_{(1)}>x_{(1)}, X_{(n)}<x_{(n)}\\ =\Phi\left(\tfrac{x_{(n)}-\mu}{\sigma}\right)^n - \left[\Phi\left(\tfrac{x_{(n)}-\mu}{\sigma}\right) -\Phi\left(\tfrac{x_{(1)}-\mu}{\sigma}\right)\right]^n$

donde es el CDF gaussiano estándar. La diferenciación con respecto a y da la función de densidad de probabilidad conjunta $\Phi(\cdot)$ $x_{(1)}$ $x_{(n)}$

F (X_{(1)}, X_{(norte)}; μ, σ) = norte (norte - 1) {[Φ (\frac{X_{(norte)} - μ}{σ}) - Φ (\frac{X_{(1)} - μ}{σ})]}^{norte - 2} \cdot ϕ (\frac{X_{(norte)} - μ}{σ}) \cdot ϕ (\frac{X_{(1)} - μ}{σ}) \cdot \frac{1}{σ^{2}}

$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 es el PDF gaussiano estándar. Tomar los términos de registro y descarte que no contienen parámetros proporciona la función de probabilidad de registro $\phi(\cdot)$

ℓ (μ, σ; x_{(1)}, x_{(n)}) = (n - 2) \log [Φ (\frac{x_{(n)} - μ}{σ}) - Φ (\frac{x_{(1)} - μ}{σ})] + \log ϕ (\frac{x_{(n)} - μ}{σ}) + \log ϕ (\frac{x_{(1)} - μ}{σ}) - 2 \log σ

$\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 se ve muy tratable pero es fácil ver que se maximiza cualquiera que sea el valor de por el ajuste $\sigma$ , es decir, el punto medio: el primer término se maximiza cuando el argumento de un CDF es negativo del argumento del otro; los términos segundo y tercero representan la probabilidad conjunta de dos variables normales independientes. $\mu=\hat\mu=\frac{x_{(n)}+x_{(1)}}{2}$

Sustituyendo en el diario de probabilidad y escribir da $\hat\mu$ $r=x_{(n)}-x_{(1)}$

ℓ (σ; x_{(1)}, x_{(n)}, \hat{μ}) = (n - 2) \log [1 - 2 Φ (\frac{- r}{2 σ})] - \frac{r^{2}}{4 σ^{2}} - 2 \log σ

$\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 tiene que ser maximizado numéricamente (por ejemplo, con el optimizede la R statpaquete) para encontrar . (Resulta que , donde es una constante que sólo depende de -tal vez a alguien más matemáticamente hábil de lo que podía demostrar por qué.) $\hat\sigma$ $\hat\sigma=k(n)\cdot r$ $k$ $n$

Las estimaciones no sirven sin una medida de precisión que lo acompañe. La información de Fisher observada puede evaluarse numéricamente (por ejemplo, con hessianel numDerivpaquete de R ) y usarse para calcular errores estándar aproximados:

I (μ) = - {\frac{\partial^{2} ℓ (μ; \hat{σ})}{(\partial μ)^{2}} |}_{μ = \hat{μ}}

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

I (σ) = - {\frac{\partial^{2} ℓ (σ; \hat{μ})}{(\partial σ)^{2}} |}_{σ = \hat{σ}}

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

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

Scortchi - Restablece a Monica
fuente

+1. Agregar la constante

a la probabilidad de log no cambiará la ubicación de su máximo, pero la convierte en una función de

, de donde el valor de

que lo maximiza es alguna función

. De manera

como usted demanda. En otras palabras, la cantidad relevante para trabajar es la relación de la desviación estándar al rango (observado), o igualmente bien su recíproco, que está estrechamente relacionado con el

2 \log (r)

$2\log(r)$

σ / r

$\sigma/r$

n

$n$

σ / r

$\sigma/r$

n \to k (n)

$n\to k(n)$

\hat{σ} = k (n) r

$\hat\sigma=k(n)r$ Rango estudiado .

whuber

@whuber: ¡Gracias! Parece obvio en retrospectiva. Incorporaré eso en la respuesta.

Scortchi - Restablece a Monica

$\mu$ $\sigma$ $R=x_{(n)} - x_{(1)}$ $99.7$

μ + 3 σ \approx x_{(n)}

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

μ - 3 σ \approx x_{(1)}

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

Restando el segundo del primero obtenemos

6 σ \approx x_{(n)} - x_{(1)} = R

$6\sigma \approx x_{(n)} - x_{(1)}= R$

\hat{σ} = \frac{1}{6} ({\bar{x}}_{(n)} - {\bar{x}}_{(1)})

$\hat \sigma = \frac 16 \Big(\bar x_{(n)} - \bar x_{(1)}\Big)$

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

Alecos Papadopoulos
fuente

That's neither a close approximation for small

n

$n$ nor an asymptotic result for large

n

$n$ .

Scortchi - Reinstate Monica

@Stortchi Bueno, no dije que es una buena estimación, pero creo que siempre es bueno tener soluciones fáciles de implementar, incluso muy difíciles, para tener una idea cuantitativa del problema en cuestión, junto con más enfoques sofisticados y eficientes como, por ejemplo, el descrito en la otra respuesta a esta pregunta.

Alecos Papadopoulos

Yo no haría carpas en "la expectativa del rango de la muestra resulta ser aproximadamente 6 veces la desviación estándar para los valores de

n

$n$ entre 200 y 1.000" Pero me estoy perdiendo algo sutil en su derivación, o ¿no funcionar igual de bien para justificar la división del rango por cualquier número.?

Scortchi - Restablecer Mónica

@Scortchi Bueno, el espíritu del enfoque es "si esperamos que casi todas las realizaciones caigan dentro de 6 sigmas, entonces es razonable esperar que las realizaciones extremas estén cerca de la frontera", eso es todo, realmente. Quizás estoy demasiado acostumbrado a operar bajo información extremadamente incompleta, y me veo obligado a decir algo cuantitativo al respecto ... :)

Alecos Papadopoulos

I could reply that even more observations would fall within

10 σ

$10 \sigma$ of the mean, giving a better estimate

\hat{σ} = \frac{R}{10}

$\hat\sigma=\frac{R}{10}$ . I shan't because it's nonsense. Any number over

1.13

$1.13$ will be a rough estimate for some value of

n

$n$ .

Scortchi - Reinstate Monica

Es sencillo obtener la función de distribución del máximo de la distribución normal (ver "P.max.norm" en el código). De él (con algunos cálculos) puede obtener la función cuantil (ver "Q.max.norm").

Usando "Q.max.norm" y "Q.min.norm" puede obtener la mediana del rango que está relacionado con N. Usando la idea presentada por Alecos Papadopoulos (en respuesta anterior) puede calcular sd.

Prueba esto:

N = 100000    # the size of the sample

# Probability function given q and N
P.max.norm <- function(q, N=1, mean=0, sd=1){
    pnorm(q,mean,sd)^N
} 
# Quantile functions given p and 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)
} 

### lets test it (takes some time)
Q.max.norm(0.5, N=N)  # The median on the maximum
Q.min.norm(0.5, N=N)  # The median on the minimum

iter = 100
median(replicate(iter, max(rnorm(N))))
median(replicate(iter, min(rnorm(N))))
# it is quite OK

### Lets try to get estimations
true_mean = -3
true_sd = 2
N = 100000

x = rnorm(N, true_mean, true_sd)  # simulation
x.vec = range(x)                  # observations

# estimation
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)

# Quite good, but only for large N
# -3  2
# -3.252606  1.981593

Vyga
fuente

Continuando con este enfoque,

E (R) = σ \int_{- \infty}^{\infty} 1 - (1 - Φ (x))^{n} - Φ (x)^{n} d x = σ d_{2} (n)

$\operatorname{E} (R) = \sigma \int_{-\infty}^{\infty} 1-(1-\Phi(x))^n -\Phi(x)^n\, \mathrm{d} x = \sigma d_2(n)$ , dónde

R

$R$ es el rango y

Φ (\cdot)

$\Phi(\cdot)$ La función de distribución acumulativa normal estándar. Puede encontrar valores tabulados de

d_{2}

$d_2$ Para pequeños

n

$n$ en la literatura de control de procesos estadísticos, evalúe numéricamente la integral o simule para su

n

$n$ .

Scortchi - Restablece a Monica