14
$\begingroup$

I know this might be a little ropey, statistically, but this is my problem.

I have a lot of range data, that is to say the minimum, maximum and sample size of a variable. For some of these data I also have a mean, but not many. I want to compare these ranges to each other to quantify the variability of each range, and also to compare the means. I have a good reason to assume that the distribution is symmetrical around the mean, and that the data will have a Gaussian distribution. For this reason I am thinking I can justify using the mid-point of the distribution as a proxy for the mean, when it is absent.

What I want to do is reconstruct a distribution for each range, and then use that to provide a standard deviation or standard error for that distribution. The only information I have is the max and min observed from a sample, and the mid-point as a proxy for the mean.

In this way I hope to be able to calculate weighted means for each group, and also to work out the coefficient of variation for each group as well, based on the range data I have and my assumptions (of a symmetrical and normal distribution).

I plan to use R to do this, so any code help would be appreciated as well.

$\endgroup$
  • 2
    $\begingroup$ I was wondering why you say you have data for minimum & maximum & maximum values; then later that you have information on only expected minimum & maximum. Which is it - observed or expected? $\endgroup$ – Scortchi Feb 11 '14 at 13:48
  • $\begingroup$ Sorry, that is my mistake. The maximum and minimum data are observed (measured from real life objects). I have amended the post. $\endgroup$ – green_thinlake Feb 12 '14 at 14:44
11
$\begingroup$

The joint cumulative distribution function for the minimum $x_{(1)}$ & maximum $x_{(n)}$ for a sample of $n$ from a Gaussian distribution with mean $\mu$ & standard deviation $\sigma$ is

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

where $\Phi(\cdot)$ is the standard Gaussian CDF. Differentiation with respect to $x_{(1)}$ & $x_{(n)}$ gives the joint probability density function

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

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

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

This doesn't look very tractable but it's easy to see that it's maximized whatever the value of $\sigma$ by setting $\mu=\hat\mu=\frac{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.

Substituting $\hat\mu$ into the log-likelihood & writing $r=x_{(n)}-x_{(1)}$ gives $$\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}$$

This expression has to be maximized numerically (e.g. with optimize from R's stat package) to find $\hat\sigma$. (It turns out that $\hat\sigma=k(n)\cdot 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(\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}$$

It would be interesting to compare the likelihood & the method-of-moments estimates for $\sigma$ 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.

$\endgroup$
  • 1
    $\begingroup$ +1. Adding the constant $2\log(r)$ to the log-likelihood will not change the location of its maximum, but converts it into a function of $\sigma/r$ and $n$, whence the value of $\sigma/r$ that maximizes it is some function $n\to k(n)$. Equivalently, $\hat\sigma=k(n)r$ as you claim. In other words, the relevant quantity to work with is the ratio of the standard deviation to the (observed) range, or equally well its reciprocal--which is closely related to the Studentized range. $\endgroup$ – whuber May 7 '14 at 19:47
  • $\begingroup$ @whuber: Thanks! Seems obvious with hindsight. I'll incorporate that into the answer. $\endgroup$ – Scortchi May 8 '14 at 7:28
1
$\begingroup$

You need to relate the range to the standard deviation/variance.Let $\mu$ be the mean, $\sigma$ 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,

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

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

Subtracting the second from the first we obtain

$$6\sigma \approx 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 $$\hat \sigma = \frac 16 \Big(\bar x_{(n)} - \bar x_{(1)}\Big)$$ 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.

$\endgroup$
  • 3
    $\begingroup$ That's neither a close approximation for small $n$ nor an asymptotic result for large $n$. $\endgroup$ – Scortchi Feb 7 '14 at 12:39
  • 1
    $\begingroup$ @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. $\endgroup$ – Alecos Papadopoulos Feb 7 '14 at 17:50
  • $\begingroup$ 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? $\endgroup$ – Scortchi Feb 7 '14 at 18:23
  • $\begingroup$ @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... :) $\endgroup$ – Alecos Papadopoulos Feb 7 '14 at 18:32
  • 4
    $\begingroup$ I could reply that even more observations would fall within $10 \sigma$ of the mean, giving a better estimate $\hat\sigma=\frac{R}{10}$. I shan't because it's nonsense. Any number over $1.13$ will be a rough estimate for some value of $n$. $\endgroup$ – Scortchi Feb 8 '14 at 0:52
1
$\begingroup$

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
$\endgroup$
  • 2
    $\begingroup$ Continuing this approach, $\operatorname{E} (R) = \sigma \int_{-\infty}^{\infty} 1-(1-\Phi(x))^n -\Phi(x)^n\, \mathrm{d} x = \sigma d_2(n)$, where $R$ is the range & $\Phi(\cdot)$ the standard normal cumulative distribution function. You can find tabulated values of $d_2$ for small $n$ in the statistical process control literature, numerically evaluate the integral, or simulate for your $n$. $\endgroup$ – Scortchi Feb 7 '14 at 14:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.