The statistical probability of state of a solution containing a reacting receptor M, a ligand A (and eventually proton H) is described by a molar partition function Z(M) = exp(-Delta G/RT) referred to M, or Z(A) or Z(H), respectively. The partition function for one class of sites can be expressed as the function of site constants k(j) and cooperativity functions gamma(j,i) = exp{b(j)(i - 1)}. Binding in a single class can be represented by a vector J(p(t)) whose elements correspond to single species. For more classes of sites, the joined probability is obtained as tenser product of single class vectors giving rise to complexation matrices M-pqr. There is one partition function for each component of the system. If the complexes are of type H(P)M(Q)A(R) there are three partition functions Z(H), ZM and ZA The relationships between partition functions and total analytical amounts T-H, T-M, T-A, respectively are given. The experimental data obtained in a potentiometric titration with electrode reversible to [H] or other free component can be reproduced as the function of site constants k(j) and cooperativity functions exp{b(j)(i - 1)} for each class j. The best values of k(j) and b(j), can be calculated following a nonlinear least squares procedure by means of a computer program that is here presented.