[R] Multiple Integrals

[Charles C. Berry]

On Fri, 28 Aug 2015, Shant Ch via R-help wrote:

> Hello all,

> For a study I want to find E|X1/3+X2/3+X3/3-X4| for variables following 
> lognormal distribution. To get the value we need to use four integrals.

So far, so good.


> This is the code which I is used
>
>       fx<-function(x){
>         dlnorm(x,meanlog=2.185,sdlog=0.562)
>       }
> U31<-integrate(function(y1) { sapply(y1, function(y1) { 
> +              integrate(function(y2){  sapply(y2, function(y2) {
> +              integrate(function(x1){  sapply(x1, function(x1) { 
> +              integrate(function(x2)
> +               abs(y1/3+y2/3+x1/3-x2)*fx(y1)*fx(y2)*fx(x1)*fx(x2),0, Inf)$value
> +              })},0, Inf)$value })},0, Inf)$value})},0,Inf)$value
> The error I received is the following:
> Error in integrate(function(y2) { :
>   maximum number of subdivisions reached
>
> I can understand the problem,

This is NOT a problem for which a numerical solution (apart from 
evaluating exp(x) and then doing some arithmetic) is required.

You are calculating the expectation of a sum of random variables. And from 
your code, these are independent random variables.

There is a well known calculus of expectations. Consult a book on 
probability theory. The expectation of a lognormal distribution is 
both well known and easy to work out. So is the expectation of a constant 
times a random variable. The expectation of a sum of random variables is 
also well known.

So write down the expectation of the lognormal. Render that expression as 
an R function.

Write down the expectation of a random variable times a constant. Render 
that expression as an R function.

Write down the expression for expectation of a sum of independent 
random variables as a function of the expectations of the random 
variables.

Lastly, write an R function that calls all of the above to yield the 
expectation of a sum of lognormal variables times fixed constants.

You do not need to use (and should not use!) the integrate() function to 
accomplish your aim.

Chuck

p.s. Nesting calls to integrate() is almost certainly a very bad approach 
to any multiple integration problem. If you ever do need to solve a 
multiple integration problem numerically, consult an expert before trying 
to write the solution as R code.


> but I am unable to figure out what can be done.. It would be great if 
> you can let me know a solution to the problem so as to find a value for 
> the integral.
>
> Shant
>
> 	[[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.

Charles C. Berry                 Dept of Family Medicine & Public Health
cberry at ucsd edu               UC San Diego / La Jolla, CA 92093-0901
http://famprevmed.ucsd.edu/faculty/cberry/