def _fast_isserli(powerlist, Sigma_x):
    """
	Get higher order mixed centralized moments of a multivariate gaussian using isserlis theorem
	http://en.wikipedia.org/wiki/Isserlis%27_theorem

	This is optimized by assuming the x_i to be independent

	:math:`E[(x-\\mu_x)^{k_x} (y-\\mu_y)^{k_y}] = E[(x-\\mu_x)^{k_x}] \\cdot E [(y-\\mu_y)^{k_y}]` => no need for isserli

	:param powerlist: list of powers of the random variables of the multivariate normal
	:param Sigma_x: The covariance matrix
	:return:
	"""
    from scipy.misc import factorial2
    if powerlist.sum() % 2 != 0:
        return 0
    for power in powerlist:
        if power % 2 != 0:
            return 0
    part = 1.0
    for i, power in enumerate(powerlist):
        part *= Sigma_x[i][i] ** (power / 2)
        part *= factorial2(power - 1, exact=True)
    return part
