Czy mogę zrekonstruować rozkład normalny na podstawie wielkości próbki oraz wartości minimalnych i maksymalnych? Mogę użyć punktu środkowego do określenia średniej

14

Wiem, że to może być trochę ryzykowne statystycznie, ale to mój problem.

Mam wiele danych zakresu, to znaczy minimalną, maksymalną i wielkość próbki zmiennej. Dla niektórych z tych danych mam również średnią, ale nie wiele. Chcę porównać te zakresy ze sobą, aby obliczyć zmienność każdego zakresu, a także porównać średnie. Mam dobry powód przypuszczać, że rozkład jest symetryczny wokół średniej i że dane będą miały rozkład Gaussa. Z tego powodu myślę, że mogę uzasadnić użycie środkowej części rozkładu jako przybliżenia średniej, gdy jest ona nieobecna.

Chcę zrekonstruować rozkład dla każdego zakresu, a następnie użyć go, aby podać odchylenie standardowe lub błąd standardowy dla tego rozkładu. Jedyne informacje, jakie mam, to maksima i min obserwowane z próbki oraz punkt środkowy jako przybliżenie średniej.

W ten sposób mam nadzieję, że będę w stanie obliczyć średnie ważone dla każdej grupy, a także opracować współczynnik zmienności dla każdej grupy, w oparciu o dane zakresu i moje założenia (rozkład symetryczny i normalny).

Planuję użyć do tego R, więc każda pomoc kodu byłaby mile widziana.

green_thinlake
źródło
2
Zastanawiałem się, dlaczego twierdzisz, że masz dane dotyczące wartości minimalnych, maksymalnych i maksymalnych; potem masz informacje o oczekiwanym minimum i maksimum. Co to jest - zaobserwowane lub oczekiwane?
Scortchi - Przywróć Monikę
Przepraszam, to mój błąd. Obserwowane są maksymalne i minimalne dane (mierzone na podstawie rzeczywistych obiektów). Poprawiłem post.
green_thinlake

Odpowiedzi:

11

Łączna funkcja skumulowanego rozkładu dla minimum i maksimum x ( n ) dla próbki n z rozkładu Gaussa ze średnią μ i odchyleniem standardowym σ wynosix(1)x(n)nμσ

F(x(1),x(n);μ,σ)=Pr(X(1)<x(1),X(n)<x(n))=Par(X(n)<x(n))-Par(X(1)>x(1),X(n)<x(n)=Φ(x(n)-μσ)n-[Φ(x(n)-μσ)-Φ(x(1)-μσ)]n

gdzie to standardowy gaussowski CDF. Zróżnicowanie względem x ( 1 ) i x ( n ) daje funkcję gęstości prawdopodobieństwa połączeniaΦ()x(1)x(n)

f(x(1),x(n);μ,σ)=n(n1)[Φ(x(n)μσ)Φ(x(1)μσ)]n2ϕ(x(n)μσ)ϕ(x(1)μσ)1σ2

gdzie to standardowy gaussowski plik PDF. Biorąc dziennik i upuszczając warunki, które nie zawierają parametrów, daje funkcję prawdopodobieństwa dziennikaϕ()

(μ,σ;x(1),x(n))=(n2)log[Φ(x(n)μσ)Φ(x(1)μσ)]+logϕ(x(n)μσ)+logϕ(x(1)μσ)2logσ

To nie wygląda bardzo łagodny, ale to łatwo zauważyć, że bez względu na to maksymalizować wartość przez ustawienie ľ = ľ = x ( n ) + x ( 1 )σ , tj. Punkt środkowy - pierwszy termin jest maksymalizowany, gdy argument jednego CDF jest ujemny od argumentu drugiego; drugi i trzeci termin reprezentują wspólne prawdopodobieństwo dwóch niezależnych zmiennych normalnych.μ=μ^=x(n)+x(1)2

Podstawiając ľ w Log-Likelihood i pisania R = x ( n ) - x ( 1 ) daje £ -l ( σ ; x ( 1 ) , x ( n ) , μ ) = ( n - 2 ) log [ 1 - 2 Φ ( - rμ^r=x(n)x(1)

(σ;x(1),x(n),μ^)=(n2)log[12Φ(r2σ)]r24σ22logσ

Wyrażenie to ma zostać zmaksymalizowane liczbowo (np optimizez R w statzestawie) w celu znalezienia σ . (Okazuje się, że Ďσ^ , gdzie k jest stałą zależności tylko od n -perhaps ktoś bardziej matematycznie zręczny niż mogłem pokazać, dlaczego).σ^=k(n)rkn

Szacunki nie mają zastosowania bez towarzyszącej mi precyzji. Obserwowane informacje Fishera można ocenić numerycznie (np. Z pakietu hessianR numDeriv) i wykorzystać do obliczenia przybliżonych błędów standardowych:

I(σ)=-2(σ; μ )

I(μ)=2(μ;σ^)(μ)2|μ=μ^
I(σ)=2(σ;μ^)(σ)2|σ=σ^

Interesujące byłoby porównanie prawdopodobieństwa i oszacowania metody momentów dla pod względem błędu (czy MLE jest spójny?), Wariancji i błędu średniej kwadratowej. Istnieje również kwestia szacowania dla tych grup, w których średnia próbki jest znana oprócz minimum i maksimum.σ

Scortchi - Przywróć Monikę
źródło
1
2log(r)σ/rnσ/rnk(n)σ^=k(n)r studentyzowane zakres .
whuber
@whuber: Dzięki! Z perspektywy czasu wydaje się to oczywiste. Włączę to do odpowiedzi.
Scortchi - Przywróć Monikę
1

You need to relate the range to the standard deviation/variance.Let μ be the mean, σ the standard deviation and R=x(n)x(1) be the range. Then for the normal distribution we have that 99.7% of probability mass lies within 3 standard deviations from the mean. This, as a practical rule means that with very high probability,

μ+3σx(n)
and

μ3σx(1)

Subtracting the second from the first we obtain

6σx(n)x(1)=R
(this, by the way is whence the "six-sigma" quality assurance methodology in industry comes). Then you can obtain an estimate for the standard deviation by
σ^=16(x¯(n)x¯(1))
where the bar denotes averages. This is when you assume that all sub-samples come from the same distribution (you wrote about having expected ranges). If each sample is a different normal, with different mean and variance, then you can use the formula for each sample, but the uncertainty / possible inaccuracy in the estimated value of the standard deviation will be much larger.

Having a value for the mean and for the standard deviation completely characterizes the normal distribution.

Alecos Papadopoulos
źródło
3
That's neither a close approximation for small n nor an asymptotic result for large n.
Scortchi - Reinstate Monica
1
@Stortchi Well, I didn't say that it is a good estimate -but I believe that it is always good to have easily implemented solutions, even very rough, in order to get a quantitative sense of the issue at hand, alongside the more sophisticated and efficient approaches like for example the one outlined in the other answer to this question.
Alecos Papadopoulos
I wouldn't carp at "the expectation of the sample range turns out to be about 6 times the standard deviation for values of n from 200 to 1000". But am I missing something subtle in your derivation, or wouldn't it work just as well to justify dividing the range by any number?
Scortchi - Reinstate Monica
@Scortchi Well, the spirit of the approach is "if we expect almost all realizations to fall within 6 sigmas, then it is reasonable to expect that the extreme realizations will be near the border" -that's all there is to it, really. Perhaps I am too used to operate under extremely incomplete information, and obliged to say something quantitative about it... :)
Alecos Papadopoulos
4
I could reply that even more observations would fall within 10σ of the mean, giving a better estimate σ^=R10. I shan't because it's nonsense. Any number over 1.13 will be a rough estimate for some value of n.
Scortchi - Reinstate Monica
1

It is straightforward to get the distribution function of the maximum of the normal distribution (see "P.max.norm" in code). From it (with some calculus) you can get the quantile function (see "Q.max.norm").

Using "Q.max.norm" and "Q.min.norm" you can get the median of the range that is related with N. Using the idea presented by Alecos Papadopoulos (in previous answer) you can calculate sd.

Try this:

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
źródło
2
Continuing this approach, E(R)=σ1(1Φ(x))nΦ(x)ndx=σd2(n), where R is the range & Φ() the standard normal cumulative distribution function. You can find tabulated values of d2 for small n in the statistical process control literature, numerically evaluate the integral, or simulate for your n.
Scortchi - Reinstate Monica