[R] Solve an equation including integral
Berend Hasselman
bhh at xs4all.nl
Fri Jan 8 21:17:05 CET 2016
> On 8 Jan 2016, at 17:47, Ravi Varadhan <ravi.varadhan at jhu.edu> wrote:
>
> There was a mistake in the previously sent function. I had hard coded the values of parameters `nu' and `exi'.
>
> Use this one:
>
> myroot <- function(t, nu, exi, alpha){
> fun <- function(t, exi, nu) (nu+t^2)^(-(nu+1)/2)*exp(((nu+1)*exi*t)/((nu+t^2)^0.5))
> res <- alpha - integrate(fun, -Inf, t, nu=nu, exi=exi)$value
> return(res)
> }
>
> uniroot(myroot, c(-2, 2), nu=2, exi=0.5, alpha=.05)
>
Given the original question, shouldn't alpha in the myroot function be replaced by (1-alpha)?
Like so
myroot <- function(t, nu, exi, alpha){
fun <- function(t, exi, nu) (nu+t^2)^(-(nu+1)/2)*exp(((nu+1)*exi*t)/((nu+t^2)^0.5))
res <- (1-alpha) - integrate(fun, -Inf, t, nu=nu, exi=exi)$value
return(res)
}
If you have to do this frequently it would be better to move the declaration of fun to outside of myroot.
Berend
> Ravi
>
>
> From: Ravi Varadhan
> Sent: Friday, January 08, 2016 11:29 AM
> To: r-help at r-project.org; 'sstoline at gmail.com' <sstoline at gmail.com>
> Subject: Re: Solve an equation including integral
>
> I think this is what you want:
>
> myroot <- function(t, nu, exi, alpha){
> fun <- function(t, exi, nu) (nu+t^2)^(-(nu+1)/2)*exp(((nu+1)*exi*t)/((nu+t^2)^0.5))
> res <- alpha - integrate(fun, -Inf, t, nu=2, exi=0.5)$value
> return(res)
> }
>
> uniroot(myroot, c(-2, 2), nu=2, exi=0.5, alpha=.05)
>
> Hope this is helpful,
> Ravi
>
> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
> Associate Professor, Department of Oncology
> Division of Biostatistics & Bionformatics
> Sidney Kimmel Comprehensive Cancer Center
> Johns Hopkins University
> 550 N. Broadway, Suite 1111-E
> Baltimore, MD 21205
> 410-502-2619
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> 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.
More information about the R-help
mailing list