I was wondering how I could have R tell me the SD (as an argument in the qnorm() built in R) for a normal distribution whose 95% limit values are already known?
As an example, I know the two 95% limit values for my normal are 158, and 168, respectively. So, in the below R code SD is shown as "x". If "y" (the answer of this simple qnorm() function) needs to be (158, 168), then can R tell me what should be x?
y <- qnorm(c(.025,.975), 163, x)
Suppose we have a Normal distribution X ~ N(mu, sigma), with unknown mean mu and unknown standard deviation sigma. And we aim to solve for mu and sigma, given two quantile equations:
Pr(X < q1) = alpha1
Pr(X < q2) = alpha2
We consider standardization: Z = (X - mu) / sigma, so that
Pr(Z < (q1 - mu) / sigma) = alpha1
Pr(Z < (q2 - mu) / sigma) = alpha2
In other words,
(q1 - mu) / sigma = qnorm(alpha1)
(q2 - mu) / sigma = qnorm(alpha2)
The RHS is explicitly known, and we define beta1 = qnorm(alpha1), beta2 = qnorm(alpha2). Now, the above simplifies to a system of 2 linear equations:
mu + beta1 * sigma = q1
mu + beta2 * sigma = q2
This system has coefficient matrix:
1 beta1
1 beta2
with determinant beta2 - beta1. The only situation for singularity is beta2 = beta1. As long as the system is non-singular, we can use solve to solve for mu and sigma.
Think about what the singularity situation means. qnorm is strictly monotone for Normal distribution. So beta1 = beta2 is as same as alpha1 = alpha2. But this can be easily avoided as it is under your specification, so in the following I will not check singularity.
Wrap up above into an estimation function:
est <- function(q, alpha) {
beta <- qnorm(alpha)
setNames(solve(cbind(1, beta), q), c("mu", "sigma"))
}
Let's have a test:
x <- est(c(158, 168), c(0.025, 0.975))
# mu sigma
#163.000000 2.551067
## verification
qnorm(c(0.025, 0.975), x[1], x[2])
# [1] 158 168
We can also do something arbitrary:
x <- est(c(1, 5), c(0.1, 0.4))
# mu sigma
#5.985590 3.890277
## verification
qnorm(c(0.1, 0.4), x[1], x[2])
# [1] 1 5