3

In my below R code, I was wondering how I could find out what is rh1 when y == 0.5?

Note that y uses atanh(rh1), which can be converted back to rh1 using tanh().

rh1 <- seq(-1, 0.1, by = 0.001)
y <- pnorm(-0.13, atanh(rh1), 0.2)
plot(rh1, y, type = "l")
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
rnorouzian
  • 7,397
  • 5
  • 27
  • 72

2 Answers2

4

Analytical solution

For a normal distribution X ~ N(mu, 0.2). We want to find mu, such that Pr (X < -0.13) = y.

Recall your previous question and my answer over there: Determine a normal distribution given its quantile information. Here we have something simpler, as there is only one unknown parameter and one piece of quantile information.

Again, we start by standardization:

   Pr {X < -0.13} = y
=> Pr { [(X - mu) / 0.2] < [(-0.13 - mu) / 0.2] } = y
=> Pr { Z < [(-0.13 - mu) / 0.2] } = y    # Z ~ N(0,1)
=> (-0.13 - mu) / 0.2 = qnorm (y)
=> mu = -0.13 - 0.2 * qnorm (y)

Now, let atanh(rh1) = mu => rh1 = tanh(mu), so in short, the analytical solution is:

tanh( -0.13 - 0.2 * qnorm (y) )

Numerical solution

It is a root finding problem. We first build the following function f, and we aim to find its root, i.e., the rh1 so that f(rh1) = 0.

f <- function (rh1, y) pnorm(-0.13, atanh(rh1), 0.2) - y

The simplest root finding method is bisection method, implemented by uniroot in R. I recommend you reading Uniroot solution in R for how we should work with it in general.

curve(f(x, 0.5), from = -1, to = 0.1); abline (h = 0, lty = 2)

We see there is a root between (-0.2, 0), so:

uniroot(f, c(-0.2, 0), y = 0.5)$root
# [1] -0.129243
Community
  • 1
  • 1
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • `uniroot` actually does something a little fancier than bisection (Brent's method: 'The function makes use of the "gold section" procedure combined with the parabolic interpolation', from http://www.netlib.org/c/brent.shar ). Not that it would matter for this problem. – Ben Bolker Jan 12 '17 at 22:54
3

Your function is monotonic so you can just create the inverse function.

rh1 <- seq(-1,.1,by=.001)
y <- pnorm(-.13,atanh(rh1),.2)

InverseFun = approxfun(y, rh1)
InverseFun(0.5)
[1] -0.1292726
G5W
  • 36,531
  • 10
  • 47
  • 80