How to compute this triple integral?

Question

I am trying to do this triple integral $$\int_{0}^{\infty }\int_{0}^{\infty }\int_{0}^{\infty }(u+w)e^{-\frac{(u+w)^2}{2}}(v+w)e^{-\frac{(v+w)^2}{2}}(u+v)e^{-\frac{(u+v)^2}{2}}e^{-(\mu +\lambda )u}e^{-2\mu v}e^{-(\mu -\lambda )w}dudvdw$$,

where $\mu$ and $\lambda$ are positive constants.

Question: how to solve this integral?

I has considered three methods to do this. The first one is using the R package R2Cuba for numerical integral.

rm(list = ls(all.names = TRUE))
library("R2Cuba")

# parameters
mu=0.5
lambda=1

integrand <- function(arg, weight) {
u <- arg[1]
v <- arg[2]
w <- arg[3]
ff= (u+w)*exp(-(u+w)^2/(2))*(v+w)*exp(-(v+w)^2/(2))*(u+v)*exp(-(u+v)^2/(2))*exp(-(mu+lambda)*u-2*mu*v-(mu-lambda)*w)
return(ff)
}    # end integrand

NDIM <-3
NCOMP <- 1

cuhre(NDIM, NCOMP, integrand, rel.tol= 1e-3, abs.tol= 1e-12,  flags= list(verbose=2, final=0))

# or by using the following function
#vegas(NDIM, NCOMP, integrand, rel.tol=1e-3,  abs.tol=1e-12)

And this code gives me the result Iteration 1: 127 integrand evaluations so far [1] 0.0557115 +- 3.83297e-006 chisq 0 (0 df) Iteration 2: 381 integrand evaluations so far [1] 0.0557115 +- 3.31385e-007 chisq 3.34279e-005 (1 df) integral: 0.05571151 (+-3.3e-07) nregions: 2; number of evaluations: 381; probability: 0.004613091

Then the second method by using Monte Carlo integration `

rm(list = ls(all.names = TRUE))
mu = 0.5
lambda = 1
N = 100000

u=rep(0,N)
v=rep(0,N)
w=rep(0,N))

sum_triple = 0 
k=0

while(k<=N) {
  u = rexp(1)
  v = rexp(1)
  w = rexp(1)
  k = k+1

  sum_triple = sum_triple + (u+w)*exp(-(u+w)^2/(2))*(v+w)*exp(-(v+w)^2/(2))*(u+v)*exp(-(u+v)^2/(2))*exp(-(mu+lambda)*u-2*mu*v-(mu-lambda)*w)/exp(-u-v-w)
}

value_of_triple = sum_triple/N
value_of_triple

`

which give me the result [1] 0.1101979

Last, I use Jacobian matrix, with $u+w=r$, $v+w=q$ and $u+v=m$. Since we need to assure that $u$, $v$ and $w$ are positive, then we have $\left | m-q \right |\leq r\leq (m+q)$. So the integral becomes

$$ \int_{0}^{\infty }\int_{0}^{\infty }\int_{\left | m-q \right |}^{m+q}re^{-\frac{r2}{2}}\cdot qe^{-\frac{q^{2}}{2}+(\lambda -\mu )q} \cdot me^{-\frac{m2}{2}-(\mu +\lambda )m}\left | -2 \right |drdqdm $$

Under change of variables, I integrate over $r$ first, which gives me $-(e^{-\frac{(m+q)^{2}}{2}}-e^{-\frac{(m-q)^{2}}{2}})$ and getting rid of the absolute value lower bound. I get $$ \frac{1}{2}\int_{0}^{\infty }\int_{0}^{\infty }(-1)(e^{-\frac{(m+q)^{2}}{2}}-e^{-\frac{(m-q)^{2}}{2}})me^{-\frac{m2}{2}-(\mu +\lambda )m}\cdot qe^{-\frac{q^{2}}{2}+(\lambda -\mu )q}dmdq,$$

then integrate over $m$, then $q$.

Am I doing anything incorrect in my codes or in the change of variables? Is there any other method to solve this?

Thank you very much.

EDIT: (1) I use MC integral to check the Jacobian transform, which gives me [1] 0.1092351. So the transform looks correct.

(2) I also find that the package R2Cuba might not work well when r.v. going to infinity.