CODEWITHSUNDEEP
r
Zhili
Zhili

Reputation: 53

Why "Uniroot" command gives the different results

I have a question on uniroot command. I cannot raise my question in a logical way, because I have no idea about why the outcomes is different each time in the second case of my example bellow. In the first case the out come is always the same for my f function:

library(mvtnorm)
f <- function(y1_upr,y2_upr = -0.05453663,target = 25e-4,df=3) {
    pmvt(upper = c(y1_upr,y2_upr),df = df) - target
}
uniroot(f,c(-10000,10000))$root

But I don't know why when I add another argument in the same function, the result changes each time (see bellow):

g <- function(y1_upr,
              y2_upr = -0.05453663,
              y3_upr = -0.06236616,
              target = 25e-4,
              df = 3) {
    pmvt(upper = c(y1_upr,y2_upr,y3_upr),df = df) - target
} 
uniroot(g,c(-10000,10000))$root

You will see that each time you apply uniroot command with g function, it gives you different result. Does anyone has an idea on that? Can I find a way to fix my result?

Upvotes: 3

Views: 702

Answers (2)

James
James

Reputation: 66844

@joran is correct in that the default GenzBretz algorithm of pmvt has a random element involved. There is an alternative algorithm, TVPACK which gives consistent results, but may be less generally applicable:

replicate(10,pmvt(upper=c(0,-0.05453663,-0.06236616),df=3))
 [1] 0.1145384 0.1145367 0.1145365 0.1145336 0.1145381 0.1145354 0.1145385
 [8] 0.1145369 0.1145384 0.1145385
replicate(10,pmvt(upper=c(0,-0.05453663,-0.06236616),df=3,algorithm=TVPACK()))
 [1] 0.1145364 0.1145364 0.1145364 0.1145364 0.1145364 0.1145364 0.1145364
 [8] 0.1145364 0.1145364 0.1145364

Upvotes: 2

joran
joran

Reputation: 173627

The algorithm probably involves some measure of randomness in order to choose a starting point. Try:

set.seed(1) 
uniroot(g,c(-10000,10000))$root

Upvotes: 3

Related Questions