  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\'evy 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.
