Understanding the thermodynamic properties of multicomponent mixtures is of critical importance in many chemical and industrial applications. Experimental measurements become progressively difficult as the number of solution components increases -- producing the need for predictive models. Problems in development of predictive models arise if the mixture has one or more components that interact through molecular complexation or association. Experimental solubilities of anthracene and pyrene dissolved in binary systems containing one or more alcohols were measured in order to address this problem. Alcohols examined in this study were: 1-propanol, 2-propanol, 1-butanol, 2-butanol, 2-methyl-1-propanol, 3-methyl-1-butanol, and 1-octanol. In binary solvent mixtures containing only a single self-associating alcoholic solvent, the alkane cosolvents studied were: n-hexane, n-heptane, n-octane, 2,2,4-trimethylpentane, cyclohexane, methylcyclohexane, tert-butylcyclohexane. Predictive solubility equations were developed using mobile order theory. This approach differs from classical solution models by representing hydrogen bonding with a probability term rather than with expressions derived from stepwise equilibria or expressions to represent hypothetical solution aggregates. Results were compared with the predicted solubilities found from using expressions developed using the Kretschmer-Wiebe and Mecke-Kempter approaches for modeling associated solutions. It was found that the mobile order approach provided reasonably accurate predictions for the solute solubilities in the systems studied. The limitations and applications for mathematical methods of representing experimental isothermal solubility data were also studied for 72 systems. Two possible descriptive forms for this mathematical representation were suggested based on the various nearly ideal binary solvent (NIBS) and modified Wilson models.