Automated microscopic analysis of stained histopathological images is degraded by the amount of color and intensity variations in data. This paper presents a novel unsupervised probabilistic approach by integrating a convolutional neural network (CNN) and the Gaussian mixture model (GMM) in a unified framework, which jointly optimizes the modeling and normalizing the color and intensity of hematoxylin- and eosin-stained (H&E) histological images. In contrast to conventional GMM-based methods that are applied only on the color distribution of data for stain color normalization, our proposal learns how to cluster the tissue structures according to their shape and appearance and simultaneously fits a multivariate GMM to the data. This approach is more robust than standard GMM in the presence of strong staining variations because fitting the GMM is conditioned on the appearance of tissue structures in the density channel of an image. Performing a gradient descent optimization in an end-to-end learning, the network learns to maximize the log-likelihood of data given estimated parameters of multivariate Gaussian distributions. Our method does not need ground truth, shape and color assumptions of image contents or manual tuning of parameters and thresholds which makes it applicable to a wide range of histopathological images. Experiments show that our proposed method outperforms the state-of-the-art algorithms in terms of achieving a higher color constancy.