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/