Kann ich eine Normalverteilung aus dem Stichprobenumfang und den Min- und Max-Werten rekonstruieren? Ich kann den Mittelpunkt verwenden, um den Mittelwert darzustellen

14

Ich weiß, dass dies statistisch gesehen vielleicht ein bisschen blöd ist, aber das ist mein Problem.

Ich habe viele Bereichsdaten, das heißt das Minimum, das Maximum und die Stichprobengröße einer Variablen. Für einige dieser Daten habe ich auch einen Mittelwert, aber nicht viele. Ich möchte diese Bereiche miteinander vergleichen, um die Variabilität jedes Bereichs zu quantifizieren und auch die Mittelwerte zu vergleichen. Ich habe einen guten Grund anzunehmen, dass die Verteilung um den Mittelwert symmetrisch ist und dass die Daten eine Gaußsche Verteilung haben werden. Aus diesem Grund denke ich, dass ich es rechtfertigen kann, den Mittelpunkt der Verteilung als Proxy für den Mittelwert zu verwenden, wenn er fehlt.

Was ich tun möchte, ist, eine Verteilung für jeden Bereich zu rekonstruieren und diese dann zu verwenden, um eine Standardabweichung oder einen Standardfehler für diese Verteilung bereitzustellen. Die einzige Information, die ich habe, ist das von einer Stichprobe beobachtete Maximum und Minimum und der Mittelpunkt als Proxy für den Mittelwert.

Auf diese Weise hoffe ich in der Lage zu sein, gewichtete Mittelwerte für jede Gruppe zu berechnen und auch den Variationskoeffizienten für jede Gruppe zu berechnen, basierend auf den Bereichsdaten, die ich habe, und meinen Annahmen (einer symmetrischen und normalen Verteilung).

Ich plane, R zu verwenden, um dies zu tun, so würde jede mögliche Code-Hilfe ebenso geschätzt.

green_thinlake
quelle
2
Ich habe mich gefragt, warum Sie sagen, dass Sie Daten für minimale und maximale und maximale Werte haben. Später haben Sie nur Informationen zum erwarteten Minimum und Maximum. Was ist es - beobachtet oder erwartet?
Scortchi - Wiedereinsetzung von Monica
Entschuldigung, das ist mein Fehler. Die maximalen und minimalen Daten werden eingehalten (gemessen an Objekten aus dem wirklichen Leben). Ich habe den Beitrag geändert.
green_thinlake

Antworten:

11

Die gemeinsame kumulative Verteilungsfunktion für das Minimum x(1) & Maximum x(n) für eine Stichprobe von n aus einer Gaußschen Verteilung mit mittlerem μ & Standardabweichung σ ist

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)=Φ(x(n)μσ)n[Φ(x(n)μσ)Φ(x(1)μσ)]n

wobei die Standard-Gaußsche CDF ist. Die Differenzierung in Bezug auf x ( 1 ) & x ( n ) ergibt die gemeinsame WahrscheinlichkeitsdichtefunktionΦ()x(1)x(n)

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

where ϕ() is the standard Gaussian PDF. Taking the log & dropping terms that don't contain parameters gives the log-likelihood function

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

This doesn't look very tractable but it's easy to see that it's maximized whatever the value of σ by setting μ=μ^=x(n)+x(1)2, i.e. the midpoint—the first term is maximized when the argument of one CDF is the negative of the argument of the other; the second & third terms represent the joint likelihood of two independent normal variates.

μ^r=x(n)x(1)

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

This expression has to be maximized numerically (e.g. with optimize from R's stat package) to find σ^. (It turns out that σ^=k(n)r, where k is a constant depending only on n—perhaps someone more mathematically adroit than I could show why.)

Estimates are no use without an accompanying measure of precision. The observed Fisher information can be evaluated numerically (e.g. with hessian from R's numDeriv package) & used to calculate approximate standard errors:

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

It would be interesting to compare the likelihood & the method-of-moments estimates for σ in terms of bias (is the MLE consistent?), variance, & mean-square error. There's also the issue of estimation for those groups where the sample mean is known in addition to the minimum & maximum.

Scortchi - Reinstate Monica
quelle
1
+1. Die Konstante hinzufügen2Log(r) to the log-likelihood ändert den Ort seines Maximums nicht, sondern konvertiert ihn in eine Funktion von σ/r und n, woher der Wert von σ/r Das maximiert es ist eine Funktion nk(n). Äquivalent dazuσ^=k(n)rwie du behauptest. Mit anderen Worten, die relevante Größe, mit der gearbeitet werden soll, ist das Verhältnis der Standardabweichung zum (beobachteten) Bereich oder ebenso sein Kehrwert - der eng mit dem studentisierten Bereich zusammenhängt .
whuber
@whuber: Thanks! Seems obvious with hindsight. I'll incorporate that into the answer.
Scortchi - Reinstate Monica
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
quelle
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
quelle
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