It is well-documented that neural responses in sensory cortices are highly variable: the same stimulus can evoke a different response on each presentation. Traditionally, this variability has been considered as noise and eliminated by using trial-averaged responses. Such averaged responses have been used almost exclusively for characterizing neural responses and mapping receptive fields with tuning curves, and accordingly, most computational theories of cortical representations have neglected or focus on unstructured Poisson-like aspects of neural variability. However, the large magnitude, characteristic spatio-temporal patterns and systematic, stimulus-dependent changes of neural variability suggest it may play a major role in sensory processing. We propose that sensory processing and learning in humans and other animals is probabilistic following the principles of Bayesian inference, and neural activity patterns represent statistical samples from a probability distribution over visual features. In this representational scheme, the set of responses at any time in a population of neurons in V1 represents a possible combination of visual features. Variability in responses arises from the dynamics that evokes population patterns with relative frequencies equal to the probability of the corresponding combination of features under the probability distribution that needs to be represented. Consequently, the average and variability of responses encode different and complementary aspects of a probability distribution: average responses encode the mean, while variability and co-variability encode higher order moments, such as variances and covariances, of the distribution. We developed a model derived from this sampling-based representational framework and showed how it can account for the most prominent hitherto unexplained features of neural variability in V1 related to changing variability and the pattern of correlations without necessarily changing mean responses. Besides providing the traditional mean responses and ting curves, the model replicates a wide range of experimental observations on systematic variations of response variability in V1 reported in the literature. These include the quenching of variability at stimulus onset measured either by membrane potential variability or by the Fano factor of spike counts, contrast-dependent and orientation-independent variability of cell responses, contrast-dependent correlations, and the close correspondence between spontaneous and evoked response distributions in the primary visual cortex. Crucially, current theories of cortical computations do not account for any of these non-trivial aspects of neural variability. The framework also makes a number of key predictions related to the time-dependent nature of the sampling-based representation. These results suggest that representations based on samples of probability distributions provide a biologically feasible new alternative to support probabilistic inferential computations of environmental features in the brain based on isy and ambiguous inputs.