The lack of a procedure to determine equilibrium thermodynamic properties of a small system interacting with a bath is frequently seen as a weakness of conventional statistical mechanics. A typical example for such a small system is a solute surrounded by an explicit solvation shell. One way to approach this problem is to enclose the small system of interest in a large bath of explicit solvent molecules, considerably larger than the system itself. The explicit inclusion of the solvent degrees of freedom is obviously limited by the available computational resources. A potential remedy to this problem is a microsolvation approach where only a few explicit solvent molecules are considered and surrounded by an implicit solvent bath. Still, the sampling of the solvent degrees of freedom is challenging with conventional grand canonical Monte Carlo methods, since no single chemical potential for the solvent molecules can be defined in the realm of small-system thermodynamics. In this work, a statistical thermodynamic model based on the grand canonical ensemble is proposed that avoids the conventional system size limitations and accurately characterizes the properties of the system of interest subject to the thermodynamic constraints of the bath. We extend an existing microsolvation approach to a generalized multibath "microstatistical" model and show that the previously derived approaches result as a limit of our model. The framework described here is universal and we validate our method numerically for a Lennard-Jones model fluid.