We carry out the first systematic investigation of inference in probabilistic models defined through submodular functions, generalizing regular binary Markov Random Fields and Determinantal Point Processes. In particular, we present a variational approach to general log-submodular and log-supermodular distributions based on sub- and supergradients. We obtain both lower and upper bounds on the log-partition function, which enables computing probability intervals for marginals, conditionals and marginal likelihoods. We also obtain fully factorized approximate posteriors, at essentially the same computational cost as ordinary submodular optimization. Our framework results in convex problems for optimizing over differentials of submodular functions, which we show how to optimally solve. Our approximation is exact at the mode (for log-supermodular distributions), and we provide bounds on the approximation quality of the log-partition function with respect to the curvature of the function. We further establish natural relations between our variational approach and the classical mean-field method from statistical physics. Exploiting additive structure in the objective leads to scalable, parallelizable message passing algorithms. Lastly, we empirically demonstrate the accuracy of our inference scheme on several submodular models arising in computer vision and network analysis. This talk is primarily based on joint work with Josip Djolonga.