4
$\begingroup$

I would like to generate a normally distributed variable of 100 data points with a known range (e.g. $10.4\text{–}16.6$) but without a known mean. To use rnorm I need mean and sd. I could estimate a mean if I assume it is the exact midpoint (mean${} = (10.4 + 16.6) \cdot 0.5 = 13.5).$ And do some similar estimations for sd (sd${} = (13.5 - 10.4) /2 = 1.55).$ And then I could use the rnorm function

rnorm (n=100, mean = 13.5, sd = 1.55)

But the resulting variable's range is outside of my initial range of $10.4\text{–}16.6.$ I would like to see if there is a better way to generate this variable and to force it to stay true to my intended range.

$\endgroup$
6
  • 9
    $\begingroup$ A normal distribution is not bounded in principle. The sample minimum and maximum aren't reproducible therefore. You need to decide on a different approach or what you're prepared to give up as a goal. $\endgroup$
    – Nick Cox
    Commented Aug 14 at 18:36
  • 3
    $\begingroup$ Some ways to interpret this request are that "known range" means either "expected range in a sample of size $n$" or possibly "a given range for a given sample size is as likely as possible." Would those be close(r) to what you're trying to ask? $\endgroup$
    – whuber
    Commented Aug 14 at 18:44
  • 1
    $\begingroup$ Thanks for this. I think the later of what you described would be ideal. $\endgroup$
    – Dendroica
    Commented Aug 14 at 19:11
  • $\begingroup$ Technically, we would want to look at the order statistics (i.e., expected minima and maxima of a Gaussian sample of 100 values) and scale accordingly ... $\endgroup$
    – Ben Bolker
    Commented Aug 14 at 19:26
  • 1
    $\begingroup$ I think more detail is needed. Is it really necessary to constrain the range to a specific value, or does it suffice to generate a sample which is "likely" to have a similar range? $\endgroup$
    – AdamO
    Commented Aug 15 at 14:44

5 Answers 5

10
$\begingroup$

I have a solution that makes theoretical sense but doesn't force your values to stay within the range; rather, it makes the expected min/max equal to your observed min and max (thus, I'd expect about 1/4 of the generated samples to lie within your range). The order statistics are the expected values of the $r$th-ranked sample in a sample of size $n$, and this question gives a good approximation.

ordermult <- function(r, n) {
    alpha <- 0.375
    mult <- qnorm((r-alpha)/(n-2*alpha+1))
    return(mult)
}

The expected maximum of a sample of 100 (the minimum is symmetric):

m <- ordermult(100, 100)    
## 2.498591

Generating values with the mean and SD scaled so that your observed mean and range match the mean and range of the reported data (making it a function because I want to re-use it):

set.seed(101)
f <- function() {
   rnorm(n=100, mean = (10.4+16.6)/2, sd = (16.6-10.4)/(2*m))
}
f()
minvals <- replicate(1000, min(f()))
summary(minvals)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  8.268  10.084  10.447  10.390  10.767  11.668 

What fraction of runs give values within your range?

within <- replicate(10000, {x <- f(); min(x) > 10.4 && max(x) < 16.6 })
mean(within)
0.2922

29% rather than the expected 25% (I guess this is from approximation of the order statistics?)

(According to another answer, using a constant of $\pi/8$ is more precise than 0.375 in the approximation, but doing that doesn't change the acceptance rate very much, so ... ??)

$\endgroup$
1
  • 4
    $\begingroup$ +1. An alternative that adheres a little more closely to the OP's comment to the question is to note that the MLE of the mean is the midrange, so all one need do is find a variance that will yield a example likely to have the desired range. $\endgroup$
    – whuber
    Commented Aug 14 at 19:57
6
$\begingroup$

There's no such thing as "normal distribution with a known range". A normal distribution has range of $\mathbb R$.

If you want a truncated normal, you can delete those points outside of the desired range after normal samples are drawn.

https://en.wikipedia.org/wiki/Truncated_normal_distribution

$\endgroup$
1
  • $\begingroup$ To clarify, by "known range", I really meant based on a published range. Thanks. $\endgroup$
    – Dendroica
    Commented Aug 14 at 19:09
6
$\begingroup$

Unlike the other answerers, I think your question is perfectly clear. You are asking about the joint distribution of an iid $N(\mu,\sigma^2)$ sample $x_1,x_2,\dots,x_n$ conditional on the order statistics $x_{(1)},x_{(n)}$. Since you don't know $\mu$ and $\sigma$ it seems reasonable to use non-informative, improper location and scale priors on $\mu$ and $\sigma$, that is, $\pi(\mu,\sigma)\propto 1/\sigma$. You can certainly sample from the marginal posterior distribution of the sample (marginalising out $\mu$ and $\sigma$) under these assumptions using MCMC.

It turns out, however, that the resulting posterior distribution is identical to the distribution generated by simply simulating iid standard normal samples followed by a suitable shift and rescaling of each sample such the sample mininum and maximum matches the observed values.

The following numerical example provide evidence that the two methods indeed are equivalent. A proof follows after the example.

par(mfrow=c(1,2))
n <- 10
min = 0
max = 1

# Method 1
logposterior <- function(par, min, max, n) {
  mu <- par[1]
  sigma <- exp(par[2])
  dnorm(min, mu, sigma, log=TRUE) +
    dnorm(max, mu, sigma, log=TRUE) + 
       (n-2)*log(pnorm(max, mu, sigma) - pnorm(min, mu, sigma))
}
chain <- MCMCpack::MCMCmetrop1R(logposterior, theta.init=c(0,0), n=n, min=min, max=max)
#> 
#> 
#> @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#> The Metropolis acceptance rate was 0.57224
#> @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
posteriorsamples <- apply(chain, 1, function(par) {
  mu <- par[1]
  sigma <- exp(par[2])
  c(0,truncnorm::rtruncnorm(n - 2, min, max, mu, sigma),1)
})
hist(posteriorsamples)

# Method 2
rescaledsamples <- replicate(1e+4, {
  x <- rnorm(n)
  x <- (x - min(x))/(max(x) - min(x))
  x <- min + (max - min)*x
})
hist(rescaledsamples)

Created on 2024-08-15 with reprex v2.1.0

The second method can be seen as a form of fiducial inference based on the idea that the data has been generated from a standard normal iid sample $z_1,z_2,\dots,z_n$ via the transformation $x_i=\mu + \sigma z_i$. Given $x_{(1)}$ and $x_{(2)}$ and knowing the joint density of the orderstatistics $z_{(1)},z_{(2)}$, the fiducial distribution of $\mu,\sigma$ is obtained by solving \begin{align} x_{(1)} &= \mu + \sigma z_{(1)} \\ x_{(n)} &= \mu + \sigma z_{(n)} \end{align} for $\mu,\sigma$ which gives \begin{align} \mu &= \frac{x_{(1)}z_{(n)}-x_{(n)}z_{(1)}}{z_{(n)}-z_{(1)}} \\ \sigma &= \frac{x_{(n)}-x_{(1)}}{z_{(n)}-z_{(1)}} \end{align}

Thinking of $z_{(1)},z_{(n)}$ as random and the observations $x_{(1)},x_{(n)}$ as fixed, the fiducial density of $\mu,\sigma$ is then given by the usual transformation formula for joint densities giving \begin{align} f(\mu,\sigma)&=f(z_{(1)},z_{(n)})\left|\frac{\partial(z_{(1)},z_{(n)}))}{\partial(\mu,\sigma)}\right| \\&\propto \phi(z_{(1)})\phi(z_{(n)})[\Phi(z_{(n)})-\Phi(z_{(1)}]^{n-2} \left|\begin{matrix} -\frac1\sigma & -\frac1{\sigma^2}(x_{(1)}-\mu) \\ -\frac1\sigma & -\frac1{\sigma^2}(x_{(n)}-\mu) \end{matrix}\right| \\&= \phi(z_{(1)})\phi(z_{(n)})(\Phi(z_{(n)})-\Phi(z_{(1)})^{n-2} \frac1{\sigma^3}(x_{(n)}-x_{(1)}) \\&\propto\underbrace{\frac1{\sigma}}_{\text{Prior}\\\pi(\mu,\sigma)} \underbrace{\frac1\sigma\phi\left(\frac{x_{(1)} - \mu}\sigma\right)\frac1\sigma\phi\left(\frac{x_{(n)} - \mu}\sigma\right)\left[\Phi(\frac{x_{(n)} - \mu}\sigma)-\Phi(\frac{x_{(1)} - \mu}\sigma)\right]^{n-2}}_{\text{Likelihood }L(\mu,\sigma)=f(x_{(1)},x_{(n)}|\mu,\sigma)}. \end{align} The fiducial density of $\mu,\sigma$ thus happen to be identical to the posterior density of $\mu,\sigma$ based on the above non-informative prior.

$\endgroup$
3
  • 2
    $\begingroup$ This is intriguing but stretches plausibility that the results are really what the OP wants. As in myth, it may be what they are unknowingly asking for. $\endgroup$
    – Nick Cox
    Commented Aug 15 at 9:46
  • $\begingroup$ @NickCox Perhaps it is a stretch, we will hopefully hear back from the OP :) $\endgroup$ Commented Aug 15 at 12:07
  • $\begingroup$ Nitpick: the expression for the likelihood must be multiplied by the indicator that $x_{(n)}$ is not less than $x_{(1)}.$ $\endgroup$
    – whuber
    Commented Aug 17 at 13:40
3
$\begingroup$

As was mentioned previously, there is no such thing as a normal distribution on a bounded range. So let's interpret your question as @whuber suggest, "a given range for a given sample size is as likely as possible.".
This amounts to reversing the normal tolerance interval (TI). A normal tolerance interval is an interval for which you have $(1-\alpha)$% confidence that no less than $p$% of the distribution will be found (at least $p$% of the values fall in the interval, $(1-\alpha)$% of the time).
The form of such a TI is $[\bar x \pm k.s]$, where $\bar x$ is the mean, s is the sd, and k is a so-called "k-factor" which depends on the sample size. You can find tables of k-factors in many places, e.g. here, here, here, or here. Now, if you look closely at these tables, you will see that they contain different values... This is because there are various derivations, approximations, etc. for computing these k factors, and while they are all "close" they are not equal (the formal, "exact" derivation for the double sided TI requires solving an equation where the unknow is inside an integral; this has to be done numerically, and iteratively; hence many published tables are approximations). You can find a R function to return such k factor here (it seems to use the "exact" method; see Details section).
In your case, your desired interval is $[13.5 \pm 3.1]$. Say you want to generate 100 datapoints, and that you want to be 95% confident that 95% of the data will fall in that interval. The k-factor for this is 2.234. So the mean to use will be 13.5, and the sd will be $\frac {3.1} {2.234} = 1.38765$. But if you generate 100 such normal values, you can expect, most of the time, a few valkues falling outside the range.
Now you can use 99% confidence, and 99% proportion (k=3.098 - remarkably close to your 3.1 TI width - coincidence?), so your sd would need to be 1.00065... If you generate 100 samples with these parameters, you should almost always get all the data inside the range.

$\endgroup$
0
$\begingroup$

Like others have said, normal distributions do not have bounds, which means if you must fit a normal, you could either just drop observations outside your range or pick a small enough SD so that samples outside the range are highly unlikely.

Since I haven't seen it mentioned, you could use a four-parameter beta distribution which is bounded, and can approximate normal if the two shape parameters are the same and >1.

$\endgroup$

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