On Dec 24, 2013, at 9:39 AM, Aur?lien Philippot wrote:
> Dear R experts,
> I am computing the following integral.
>
> [image: \int_{1,100} \frac{1}{x+Max(x-50,0)} g(x)dx], where g is the
> density of the standard normal, and [1,100] is the domain.
>
> 1) I use the following code which works fine:
> integrand1<- function(x){1/x*dnorm(x)}
> integrand2<- function(x){1/(2*x-50)*dnorm(x)}
>
> res1<- integrate(integrand1, lower=1, upper=50)$value+
> integrate(integrand2, lower=50, upper=100)$value
> res1
> 0.1116
>
> In other words, I split the max function depending on the value of x in the
> domain.
>
> 2) Alternatively, I can also compute it by vectorizing the max function
> integrand<- function (x) ifelse(x<50, 1/x*dnorm(x) ,
1/(2*x-50)*dnorm(x))
> res4<- integrate(integrand, lower=1, upper=100)$value
> res4
> 0.1116
>
> 3) However, in both cases, the syntax is a little bit heavy and not very
> convenient if I want to add more integrals, all of them with a max in the
> integrand. Is there a way to have a more concise syntax?
> Using max or pmax directly in the definition of the integrand does not
work.
> For example:
> integrand<- function(x) {1/(x+max(x-50,0)*dnorm(x))}
Is boolean math more to your liking?
integr_both<- function(x){
(x < 50)*(1/x*dnorm(x))+
(x>= 50 & x<100)*( 1/(2*x-50)*dnorm(x) ) }
res_b<- integrate(integr_both, lower=1, upper=100)$value
> res_b
[1] 0.1116587
--
David
>
> res<- integrate(integrand, lower=1, upper=100)$value
> res
> 4.60517
>
> Thank you for any suggestion, and merry Christmas!
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
David Winsemius
Alameda, CA, USA