There are several methods. [Here](https://arxiv.org/pdf/1908.06736.pdf) is a recent one. I implemented it in Python for $n=3$ (and for a polynomial, not only a monomial):

```python
from math import factorial
from sympy import Poly
from sympy.abc import x, y, z
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    return coef / np.prod(list(range(4, 4+j)))

def integral(P, v1, v2, v3, v4):
    X = Poly(x, x, y, z, domain="RR")
    Y = Poly(y, x, y, z, domain="RR")
    Z = Poly(z, x, y, z, domain="RR")
    B = np.column_stack((v1-v4, v2-v4, v3-v4))
    newx = B[0,0]*X + B[0,1]*Y + B[0,2]*Z + v4[0]
    newy = B[1,0]*X + B[1,1]*Y + B[1,2]*Z + v4[1]
    newz = B[2,0]*X + B[2,1]*Y + B[2,2]*Z + v4[2]
    let = {x: newx.as_expr(), y: newy.as_expr(), z: newz.as_expr()}
    Q = P.subs(let, simultaneous=True).as_expr().as_poly(x, y, z)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / 6 * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])

# polynomial to integrate    
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, v1, v2, v3, v4)
print(val)
``` 

___

# EDIT

For an arbitrary $n$:

```python
from math import factorial
from sympy import Poly
import numpy as np

def term(Q, monom):
    coef = Q.coeff_monomial(monom)
    powers = list(monom)
    j = sum(powers)
    if j == 0:
        return coef
    coef = coef * np.prod(list(map(factorial, powers)))
    n = len(monom)
    return coef / np.prod(list(range(n+1, n+j+1)))

def integral(P, S):
    gens = P.gens
    n = len(gens)
    dico = {}
    v = S[n,:]
    columns = []
    for i in range(n):
        columns.append(S[i,:] - v)    
    B = np.column_stack(tuple(columns))
    for i in range(n):
        newvar = v4[i]
        for j in range(n):
            newvar = newvar + B[i,j]*Poly(gens[j], gens, domain="RR")
        dico[gens[i]] = newvar.as_expr()
    Q = P.subs(dico, simultaneous=True).as_expr().as_poly(gens)
    monoms = Q.monoms()
    s = 0.0
    for monom in monoms:
        s = s + term(Q, monom)
    return np.abs(np.linalg.det(B)) / factorial(n) * s

###############################################################################
# tetrahedron vertices
v1 = np.array([1.0, 1.0, 1.0]) 
v2 = np.array([2.0, 2.0, 3.0]) 
v3 = np.array([3.0, 4.0, 5.0]) 
v4 = np.array([3.0, 2.0, 1.0])
S = np.array([v1, v2, v3, v4])

# polynomial to integrate
from sympy.abc import x, y, z
P = Poly(x**4 + y + 2*x*y**2 - 3*z, x, y, z, domain = "RR")

# integral
val = integral(P, S)
print(val)
```

___

**Julia** implementation:

```julia
using TypedPolynomials
using LinearAlgebra

function integratePolynomialOnSimplex(P, S)
    gens = variables(P)
    n = length(gens)
    v = S[end]    
    B = Array{Float64}(undef, n, 0)
    for i in 1:n
        B = hcat(B, S[i] - v)
    end
    Q = P(gens => v + B * vec(gens))
    s = 0.0
    for t in terms(Q)
        coef = TypedPolynomials.coefficient(t)
        powers = TypedPolynomials.exponents(t)
        j = sum(powers)
        if j == 0
            s = s + coef
            continue
        end
        coef = coef * prod(factorial.(powers))
        s = s + coef / prod((n+1):(n+j))
    end
    return abs(LinearAlgebra.det(B)) / factorial(n) * s
end


#### --- EXAMPLE --- ####

# define the polynomial to be integrated
@polyvar x y z
P = x^4 + y + 2*x*y^2 - 3*z

#= be careful
if your polynomial does not involve one of the variables, 
e.g. P(x,y,z) = x⁴ + 2xy², it must be defined as a 
polynomial in x, y, and z; you can do:
@polyvar x y z
P = x^4 + 2*x*y^2 + 0.0*z
then check that variables(P) returns (x, y, z)
 =#

# simplex vertices
v1 = [1.0, 1.0, 1.0] 
v2 = [2.0, 2.0, 3.0] 
v3 = [3.0, 4.0, 5.0] 
v4 = [3.0, 2.0, 1.0]
# simplex
S = [v1, v2, v3, v4]

# integral
integratePolynomialOnSimplex(P, S)
```