Perhaps you could transform the surface integral to a volume integral using something like Gauss' law and then use NIntegrate, which works for an inequality constraint inside Boole.

No, it takes about 10 seconds before giving the DiscretizeRegion error. I'm running version 10.2. Could that be the point?

If it takes too long, it might work better to use NIntegrate with the MonteCarlo (or any of its daughters) method of integration.

Tomorrow I'll install 10.3 and try.

@Frank. Yes, I don't know why it gives that error. The division by 0 occurs only at the origin (I consider only positive values of $x$, $y$ and $z$), but then it is 0/0 and the limit is 1 for the whole region. The origin is part of the surface only when $a==1$.

If you use RegionPlot3D and change the $==$ to $\leq$, it works. The region looks like this (for $a=0.2$ for instance):

I just assumed that the surface ( $==$) is just the "skin" of this region (the orange surface). Isn't it the case?

BTW: the region is a line ( $x=y=z$) for $a=0$ and the whole first octant when $a=1$.

Thank you Frank. I think you might be right.

I rephrase your simple case as follows:

In[84]:= 
reg = ImplicitRegion[x^2 + y^2 + z^2 == 1, {x, y, z}]

    Out[84]= ImplicitRegion[x^2 + y^2 + z^2 == 1, {x, y, z}]

    In[85]:= NIntegrate[1, {x, y, z} \[Element] reg]

    Out[85]= 12.5664

And it seems to work ( $4\pi$ is the answer). But if I do the same with my integral (for instance for $a=0.5$):

reg = ImplicitRegion[Sqrt[1/3 ((x - y)^2/(x + y)^2 + (x - z)^2/(x + z)^2 + (y - z)^2/(y + z)^2)] == .5 && x >= 0 && y >= 0 && z >= 0, {x, y, z}]

and then

NIntegrate[Sqrt[
 12 ((x - y)^2/(x + y)^2 + (x - z)^2/(x + z)^2 + (y - z)^2/(y + 
      z)^2)]/(E^(x + y + z)
   Sqrt[(-((4 x (x - z))/(x + z)^3) - (4 y (y - z))/(y + z)^3)^2 + ((
     4 z (y - z))/(y + z)^3 - (4 x (x - y))/(x + y)^3)^2 + ((
     4 y (x - y))/(x + y)^3 + (4 z (x - z))/(x + z)^3)^2]), {x, y, 
   z} \[Element] reg]

Then I get the error: "DiscretizeRegion was unable to discretize the region \ ImplicitRegion[<<2>>]. "

Any further ideas?

Hi Key, thanks for your interest. Yes, it exists. It is actually a density function for random variable $a$. The $==$ is in the Boole command, representing the region of integration, a surface defined by $\sqrt{\frac{1}{3} ((\frac{x-y}{x+y})^2+(\frac{y-z}{y+z})^2+(\frac{x-z}{x+z})^2)}=a$.

Numerical integration on an infinite domain is possible with MonteCarlo techniques.