Puis-je reconstruire une distribution normale à partir de la taille de l'échantillon et des valeurs min et max? Je peux utiliser le point médian pour représenter la moyenne

Je sais que cela pourrait être un peu compliqué, statistiquement, mais c'est mon problème.

J'ai beaucoup de données de plage, c'est-à-dire la taille minimum, maximum et échantillon d'une variable. Pour certaines de ces données, j'ai également une moyenne, mais pas beaucoup. Je veux comparer ces gammes entre elles pour quantifier la variabilité de chaque gamme, et aussi comparer les moyennes. J'ai une bonne raison de supposer que la distribution est symétrique autour de la moyenne et que les données auront une distribution gaussienne. Pour cette raison, je pense pouvoir justifier l'utilisation du point médian de la distribution comme proxy de la moyenne, lorsqu'elle est absente.

Ce que je veux faire, c'est reconstruire une distribution pour chaque plage, puis l'utiliser pour fournir un écart-type ou une erreur standard pour cette distribution. Les seules informations dont je dispose sont les valeurs max et min observées à partir d'un échantillon et le point médian comme indicateur de la moyenne.

De cette façon, j'espère être en mesure de calculer des moyennes pondérées pour chaque groupe, ainsi que de calculer également le coefficient de variation pour chaque groupe, sur la base des données de plage que j'ai et de mes hypothèses (d'une distribution symétrique et normale).

Je prévois d'utiliser R pour ce faire, donc toute aide de code serait également appréciée.

r normal-distribution estimation missing-data order-statistics green_thinlake
la source

Je me demandais pourquoi vous dites que vous avez des données pour les valeurs minimum & maximum & maximum; puis plus tard, vous avez des informations sur seulement le minimum et le maximum attendus. Quel est-il - observé ou attendu?

Scortchi - Réintégrer Monica

Désolé, c'est mon erreur. Les données maximales et minimales sont observées (mesurées à partir d'objets réels). J'ai modifié le poste.

green_thinlake

Réponses:

La fonction de distribution cumulative conjointe pour le minimum $x_{(1)}$ et le maximum $x_{(n)}$ pour un échantillon de $n$ d'une distribution gaussienne avec la moyenne $\mu$ et l'écart type $\sigma$ est

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$

où $\Phi(\cdot)$ est le CDF gaussien standard. La différenciation par rapport à $x_{(1)}$ et $x_{(n)}$ donne la fonction de densité de probabilité conjointe

f (x_{(1)}, x_{(n)}; μ, σ) = n (n - 1) {[Φ (\frac{x_{(n)} - μ}{σ}) - Φ (\frac{x_{(1)} - μ}{σ})]}^{n - 2} \cdot ϕ (\frac{x_{(n)} - μ}{σ}) \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}$

où $\phi(\cdot)$ est le PDF gaussien standard. Prendre les termes log & drop qui ne contiennent pas de paramètres donne la fonction log-vraisemblance

ℓ (μ, σ; 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$

Cela ne semble pas très traitable mais il est facile de voir qu'il est maximisée quelle que soit la valeur de par le réglage $\sigma$ , c'est-à-dire le point médian - le premier terme est maximisé lorsque l'argument d'un CDF est le négatif de l'argument de l'autre; les deuxième et troisième termes représentent la probabilité conjointe de deux variables normales indépendantes. $\mu=\hat\mu=\frac{x_{(n)}+x_{(1)}}{2}$

Substituant dans la log-vraisemblance et écriture donne $\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}$

Cette expression doit être numériquement maximisé (par exemple avec optimizede R de statpaquet) pour trouver . (Il se trouve que , où est une constante ne dépendant que de -perhaps quelqu'un mathématiquement plus adroite que je ne pouvais montrer pourquoi.) $\hat\sigma$ $\hat\sigma=k(n)\cdot r$ $k$ $n$

Les estimations ne sont pas utiles sans une mesure de précision qui l'accompagne. Les informations de Fisher observées peuvent être évaluées numériquement (par exemple avec hessianle numDerivpackage de R ) et utilisées pour calculer des erreurs standard approximatives:

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}$

Il serait intéressant de comparer la probabilité et les estimations de la méthode des moments pour en termes de biais (la MLE est-elle cohérente?), De variance et d'erreur quadratique moyenne. Il y a aussi la question de l'estimation pour les groupes où la moyenne de l'échantillon est connue en plus du minimum et du maximum. $\sigma$

Scortchi - Réintégrer Monica
la source

+1. Ajouter la constante

à la log-vraisemblance ne changera pas l'emplacement de son maximum, mais la convertit en une fonction de

, d'où la valeur de

qui la maximise est une fonction

comme vous le prétendez. En d'autres termes, la quantité pertinente à utiliser est le rapport de l'écart-type à la plage (observée), ou tout aussi bien sa réciproque - qui est étroitement liée à la

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$ Gamme étudiée .

whuber

@whuber: Merci! Semble évident avec le recul. Je vais incorporer cela dans la réponse.

Scortchi - Réintégrer Monica

Vous devez relier la plage à l'écart-type / variance. Soit la moyenne, l'écart-type et la plage. Ensuite, pour la distribution normale, nous avons que % de la masse de probabilité se situe à moins de 3 écarts-types de la moyenne. Cela signifie, en règle générale, qu'avec une probabilité très élevée, $\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)}$

En soustrayant le deuxième du premier, nous obtenons

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)$

Le fait d'avoir une valeur pour la moyenne et pour l'écart type caractérise complètement la distribution normale.

Alecos Papadopoulos
la source

n

$n$

n

$n$

@Stortchi Eh bien, je n'ai pas dit que c'était une bonne estimation, mais je pense qu'il est toujours bon d'avoir des solutions faciles à mettre en œuvre, même très approximatives, afin d'avoir une idée quantitative du problème en question, à côté des autres. des approches sophistiquées et efficaces comme par exemple celle décrite dans l'autre réponse à cette question.

Alecos Papadopoulos

Je ne me moquerais pas de "l'attente de la plage d'échantillonnage se révèle être environ 6 fois l'écart-type pour les valeurs de

n

$n$ de 200 à 1000 ". Mais est-ce que je manque quelque chose de subtil dans votre dérivation, ou cela ne fonctionnerait-il pas aussi bien pour justifier la division de la plage par un nombre quelconque?

Scortchi - Reinstate Monica

@Scortchi Eh bien, l'esprit de l'approche est "si nous nous attendons à ce que presque toutes les réalisations se situent dans les 6 sigmas, alors il est raisonnable de s'attendre à ce que les réalisations extrêmes soient près de la frontière" - c'est tout ce qu'il y a vraiment. Peut-être que je suis trop habitué à opérer sous des informations extrêmement incomplètes, et obligé de dire quelque chose de quantitatif à ce sujet ... :)

Alecos Papadopoulos

Je pourrais répondre que davantage d’observations

10 σ

$10 \sigma$ de la moyenne, donnant une meilleure estimation

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

$\hat\sigma=\frac{R}{10}$ . Je ne le serai pas parce que c'est un non-sens. Tout nombre supérieur

1.13

$1.13$ sera une estimation approximative d'une certaine valeur de

n

$n$ .

Scortchi - Réintégrer Monica

Il est simple d'obtenir la fonction de distribution du maximum de la distribution normale (voir "P.max.norm" dans le code). De là (avec un peu de calcul), vous pouvez obtenir la fonction quantile (voir "Q.max.norm").

En utilisant "Q.max.norm" et "Q.min.norm" vous pouvez obtenir la médiane de la plage qui est liée à N. En utilisant l'idée présentée par Alecos Papadopoulos (dans la réponse précédente), vous pouvez calculer sd.

Essaye ça:

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
la source

Poursuivant cette approche,

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)$ , where

R

$R$ is the range &

Φ (\cdot)

$\Phi(\cdot)$ the standard normal cumulative distribution function. You can find tabulated values of

d_{2}

$d_2$ for small

n

$n$ in the statistical process control literature, numerically evaluate the integral, or simulate for your

n

$n$ .

Scortchi - Reinstate Monica