[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