This paper proposes a Sum of Squares (SOS) optimization technique for using multivariate data to estimate the probability density function of a non-Gaussian generating process. The class of distributions over which we optimize, result from using a polynomial map to lift the data into a higher-dimensional space, solving for an optimal Gaussian fit in this space, and then projecting a polynomial slice of the resulting joint density into physical space. The resulting distribution, to be called Sliced Normal, is able to characterize multimodal responses and strong parameter dependencies. We investigate several formulations of the problem, first maximizing a log-likelihood function, then a worst-case log-likelihood function, and finally using a heuristic to increase sparsity within the maximum log-likelihood formulation - thereby identifying independent subsets of the random variables. Using the optimal density functions in each scenario, we then propose semi-algebraic sets representing confidence regions, or 'safe sets,' for future data. Finally, we show numerically that these 'safe sets' are reliable and, hence, can be used for system identification, fault detection, robustness analysis, and robust control design.