Analog behavioral models are widely used to reduce the complexity in hierarchical analog circuit design and verification. In the presence of process variations and atomic-level fluctuations, however, these models have to be extended to take variability into account. In this paper, we present a probabilistic solution that treats the behavioral model coefficients as multidimensional random variables and supports non-Gaussian as well as correlated parameters. A voltage divider and a bandgap voltage reference demonstrate the capabilities of our modeling approach in terms of accuracy and efficiency.