CODEWITHSUNDEEP
r
TourEiffel
TourEiffel

Reputation: 4424

R Solve equation - Error in uniroot(f, x = 2.326348, lower = 0.1, upper = 1e+08)

I am trying to solve in R the following equation :

2.326348 = Normal inverse function(y)

Normal Inverse Function is :

sqrt(-2*log(sqrt(2*pi)*y))

or more simply : qnorm(y)

So in order to solve the equation I used the function below :

f <- function(x,y)  (sqrt(-2*log(sqrt(2*pi)*y))-x)
uniroot(f, x=2.326348, lower=0.1, upper=100000000)$root

or

f <- function(x,y)  (qnorm(y)-x)
    uniroot(f, x=2.326348, lower=0.1, upper=100000000)$root

but this does return :

Error in uniroot(f, x = 2.326348, lower = 0.1, upper = 1e+08) : 
  f.upper = f(upper) is NA

What should I do in order to solve this equation ? I already saw this but it seems that its not the same issue.

uniroot documentation

Upvotes: 0

Views: 285

Answers (2)

G. Grothendieck
G. Grothendieck

Reputation: 269852

The question refers to both f and g below as if they were the same but they are not. The solution to f(x, y) = 0 for fixed x is y = dnorm(x) and the solution to g(x, y) = 0 for fixed x is y = pnorm(x). The expression to the left hand side of -x in the body of f is the inverse of the gaussian probability density function (pdf) while that in g is the inverse of the gaussian cumulative probability distribution function (cdf).

f <- function(x, y) (sqrt(-2*log(sqrt(2*pi)*y))-x)
g <- function(x, y) qnorm(y) - x

x <- 2.326348

f(x, dnorm(x))
## [1] 0

g(x, pnorm(x))
## [1] 4.440892e-16

1. f

First we work with f above.

Square both sides and then it will behave better in uniroot. Also use 0 and 1 as the bounds.

x <- 2.326348
f2 <- function(x,y) -2*log(sqrt(2*pi)*y)-x^2
result <- uniroot(f2, x = x, lower = 0, upper = 1)
result
## $root
## [1] 0.02666295
##
## $f.root
## [1] -0.000811161
##
## $iter
## [1] 9
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 6.103516e-05

Compare the root above with

dnorm(x)
## [1] 0.02665213

2. g

Now we work with g above.

result2 <- uniroot(g, x = x, lower = 0, upper = 1)
result2
## $root
## [1] 0.99001
##
## $f.root
## [1] 0.0003750975
##
## $iter
## [1] 11
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 6.103516e-05

Compare the root above with pnorm(x)

pnorm(x)
## [1] 0.99

Upvotes: 2

Waldi
Waldi

Reputation: 41240

The definition interval for qnorm is [0,1] and you're getting an error because you try to find a solution out of this interval.
If you limit the search to [0,1] it works :

y = 2.326348
uniroot(function(x) qnorm(x)-y,lower=0, upper = 1)
#> $root
#> [1] 0.99001
#> 
#> $f.root
#> [1] 0.0003750975
#> 
#> $iter
#> [1] 11
#> 
#> $init.it
#> [1] NA
#> 
#> $estim.prec
#> [1] 6.103516e-05

Upvotes: 1

Related Questions