Here is one solution.
import numpy as np
from scipy.integrate import tplquad
# Torus radius, cross sectional radius
R, r = 4, 1 # R, r
# Integral limits defining the torus in cylindrical coordinates
# theta limits
a, b = 0, 2 * np.pi
# r limits
gfun, hfun = lambda theta: R-r, lambda theta: R+r
# z limits
qfun = lambda theta, rho: 0
rfun = lambda theta, rho: np.sqrt(r**2 - (rho-R)**2)
lims = (a, b, gfun, hfun, qfun, rfun)
# Torus volume
V, _ = tplquad(lambda z, rho, theta: rho, *lims)
V *= 2
print('V = {} (exact: {})'.format(V, 2 * np.pi**2 * R * r**2))
# Integrand for Ix
f = lambda z, rho, theta: ((rho * np.sin(theta))**2 + z**2) * rho
Ix, _ = tplquad(f, *lims)
Ix *= 2 / V
print('Ix = {} (exact: {})'.format(Ix, 0.5 * R**2 + 5/8 * r**2))
Iz, _ = tplquad(lambda z, rho, theta: rho**3, *lims)
Iz *= 2 / V
print('Iz = {} (exact: {})'.format(Iz, R**2 + 0.75*r**2))
Output:
V = 78.95683520871499 (exact: 78.95683520871486)
Ix = 8.624999999999678 (exact: 8.625)
Iz = 16.749999999999808 (exact: 16.75)