A new family of polynomials, called cumulant polynomial sequence, and its extensions to the multivariate case is introduced relied on a purely symbolic combinatorial method. The coefficients of these polynomials are cumulants, but depending on what is plugged in the indeterminates, either sequences of moments either sequences of cumulants can be recovered. The main tool is a formal generalization of random sums, also with a multivariate random index and not necessarily integer-valued. Applications are given within parameter estimations, Lévy processes and random matrices and, more generally, problems involving multivariate functions. The connection between exponential models and multivariable Sheffer polynomial sequences offers a different viewpoint in characterizing these models. Some open problems end the paper.