Configuration entropy is believed to stabilize disordered solid solution phases in multicomponent systems at elevated temperatures over intermetallic compounds by lowering the Gibbs free energy. Traditionally the increment of configuration entropy with temperature was computed by time consuming thermodynamic integration methods. In this work, a new formalismis developed to predict configuration entropy as function of temperature frommulti-body cluster probability in a system with arbitrary number of components, K and arbitrary average composition. The formalism uses the principles behind the Cluster Variation Method (CVM) but extends the treatment to arbitrary cluster interaction employed in the Cluster Expansion (CE) technique. The multi-body probabilities are worked out by explicit inversion and direct product of a matrix formulation from point functions in the clusters obtained from symmetry independent correlation functions. The matrix quantities are determined from semi canonical Monte Carlo simulations with Effective Cluster Interactions (ECIs) derived from Density Functional Theory (DFT) calculations. The formalism is applied to analyze the four-body cluster probabilities for the quaternary system Cr-Fe-Mn-Ni. It is shown that for two specific compositions (Cr25Fe25Mn25Ni25 and Cr18Fe27Mn27Ni28) the high value of probabilities of Cr-Fe-Fe-Fe and Mn-Mn-Ni-Ni are strongly correlated with the presence of the ordering phases L12-CrFe3 and L10-MnNi respectively and that are also in excellent agreement with ab initio predictions of these ground state structures. The general formalism is used to investigate configuration entropy as a function of temperature and for 285 different alloy compositions. It is found that our matrix formulation of configuration entropy is in agreement with the result obtained by the thermodynamic integration method. At the high temperature limit, our theoretical prediction for configuration entropy agrees well with random configuration entropy for all compositions of the alloy.