R help - Apr 2019 - Vectorised uniroot function

A vectorised uniroot function would be useful for function inversion,
e.g. for quantile functions and random number generation. To address
this, I have implemented rstpm2::vuniroot that adapts the C function
R_zeroin2 for Brent's method for a vectorised objective. The function
currently uses Rcpp, but could be re-implemented using the C API.

As an example, we can now rapidly sample from a proportional hazards
mixture Weibull distribution:

pweibullMixturePH <- function(q, p1, RR, shape1, shape2, scale1=1, scale2=1)
    1 - (p1*pweibull(q, shape1, scale1, lower.tail=FALSE) +
         (1-p1)*pweibull(q, shape2, scale2, lower.tail=FALSE))^RR
rfun <- function(pfun) function(n, ..., lower=1e-5, upper=1e6) {
    u <- runif(n)
    objective <- function(q) pfun(q, ...) - u
    rstpm2::vuniroot(objective, lower=rep(lower,length=n),
                     upper=rep(upper,length=n))$root
}
rweibullMixturePH <- rfun(pweibullMixturePH)
set.seed(12345)
y <- rweibullMixturePH(n=1e4,p1=0.5,RR=2,shape1=1.5,shape2=0.5)

Has anyone previously developed a similar vectorised uniroot function?
Finally, would this be a useful addition to core R?

-- Mark


N?r du skickar e-post till Karolinska Institutet (KI) inneb?r detta att KI
kommer att behandla dina personuppgifter. H?r finns information om hur KI
behandlar
personuppgifter<https://ki.se/medarbetare/integritetsskyddspolicy>.


Sending email to Karolinska Institutet (KI) will result in KI processing your
personal data. You can read more about KI?s processing of personal data
here<https://ki.se/en/staff/data-protection-policy>.

I've had to do something similar for some of my engineering calculations.  I
would welcome something like this.  It would make the language more amenable for
engineering usage.

Thank you kindly!

Shawn Way, PE


-----Original Message-----
From: R-help <r-help-bounces at r-project.org> On Behalf Of Mark Clements
Sent: Tuesday, April 23, 2019 7:12 AM
To: r-help at r-project.org
Subject: [R] Vectorised uniroot function

** External Email **: This email originated from outside of the organization. Do
not click links or open attachments unless you recognize the sender and know the
content is safe.



A vectorised uniroot function would be useful for function inversion, e.g. for
quantile functions and random number generation. To address this, I have
implemented rstpm2::vuniroot that adapts the C function
R_zeroin2 for Brent's method for a vectorised objective. The function
currently uses Rcpp, but could be re-implemented using the C API.

As an example, we can now rapidly sample from a proportional hazards mixture
Weibull distribution:

pweibullMixturePH <- function(q, p1, RR, shape1, shape2, scale1=1, scale2=1)
    1 - (p1*pweibull(q, shape1, scale1, lower.tail=FALSE) +
         (1-p1)*pweibull(q, shape2, scale2, lower.tail=FALSE))^RR rfun <-
function(pfun) function(n, ..., lower=1e-5, upper=1e6) {
    u <- runif(n)
    objective <- function(q) pfun(q, ...) - u
    rstpm2::vuniroot(objective, lower=rep(lower,length=n),
                     upper=rep(upper,length=n))$root } rweibullMixturePH <-
rfun(pweibullMixturePH)
set.seed(12345)
y <- rweibullMixturePH(n=1e4,p1=0.5,RR=2,shape1=1.5,shape2=0.5)

Has anyone previously developed a similar vectorised uniroot function?
Finally, would this be a useful addition to core R?

-- Mark


N?r du skickar e-post till Karolinska Institutet (KI) inneb?r detta att KI
kommer att behandla dina personuppgifter. H?r finns information om hur KI
behandlar
personuppgifter<https://ki.se/medarbetare/integritetsskyddspolicy>.


Sending email to Karolinska Institutet (KI) will result in KI processing your
personal data. You can read more about KI?s processing of personal data
here<https://ki.se/en/staff/data-protection-policy>.

______________________________________________
R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.