We develop a new approach for estimating the expected values of nonlinear functions applied to multivariate random variables with arbitrary distributions. Rather than assuming a particular distribution, we assume that we are only given the first four moments of the distribution. The goal is to efficiently represent the distribution using a small number of quadrature nodes which are called $\sigma$-points. What we mean by this is choosing nodes and weights in order to match the specified moments of the distribution. The classical scaled unscented transform (SUT) matches the mean and covariance of a distribution. In this paper, we introduce the higher order unscented transform (HOUT), which also matches any given skewness and kurtosis tensors. It turns out that the key to matching the higher moments is the tensor CANDECOMP/PARAFAC (CP) decomposition. While the minimal CP decomposition is NP-complete, we present a practical algorithm for computing a nonminimal CP decomposition and prove convergence in linear time. We then show how to combine the CP decompositions of the moments in order to form the $\sigma$-points and weights of the HOUT. By passing the $\sigma$-points through a nonlinear function and applying our quadrature rule we can estimate the moments of the output distribution. We prove that the HOUT is exact on arbitrary polynomials up to fourth order and derive error bounds in terms of the regularity of the function and the decay of the probability. Finally, we numerically compare the HOUT to the SUT on nonlinear functions applied to non-Gaussian random variables including an application to forecasting and uncertainty quantification for chaotic dynamics.