Triple Integral help

Question

I cannot figure out why I'm getting a negative value running NIntegrate to calculate my triple integral. Below is my code:

S = 40
p1 = 0.37
p2 = 0.43
p3 = .05
UND = 1 - (p1 + p2 + p3)
b1 = S*p1
b2 = S*p2
b3 = S*p3
b4 = S*UND
x1 = 48
x2 = 47
x3 = 4
x4 = 1

NIntegrate[(x^(x1 + b1 - 1))*(y^(x2 + b2 - 1))*(z^(x3 + b3 - 1))*((1 -
   x - y - z)^(x4 + b4 - 1)), {z, 0, 1}, {x, z, (1 - z)/2}, {y, 0, x}]

I want to triple integrate the function:

f(x,y,z) = (x^(x1 + b1 - 1))*(y^(x2 + b2 - 1))*(z^(x3 + b3 - 1))*((1 -x - y - z)^(x4 + b4 - 1))

I want to do the following:

S S S f(x,y,z) dy dx dz

where: y goes from 0 to x

x goes from z to (1-z)/2

z goes from 0 to 1

Can anyone tell me what I'm doing wrong please?

Answer 1

As @KraZug points out:

Plot[{Callout[z, "z", {.9, Above}], 
  Callout[(1 - z)/2, "(1-z) / 2", {.1, Above}]}, {z, 0, 1}, 
 AxesLabel -> {"Z", None}]

In[2]:= Reduce[(1-z)/2>z]
Out[2]= z<1/3

So for z>1/3, the second integration range reverses.