Recent advances in machine-learning-augmented classical density functional theory (cDFT) have demonstrated that data-driven corrections to physically motivated free-energy functionals can substantially improve predictive accuracy for interfacial thermodynamics while preserving variational structure and thermodynamic structure. In particular, a recently introduced neural-network-augmented cDFT framework showed that compact, physics-informed corrections to the excess Helmholtz free energy can be calibrated against molecular dynamics data using adjoint-based optimization, enabling accurate prediction outside the training range for density profiles, fluid properties, and contact angles across scales. In this work, we extend this framework to a fully Bayesian uncertainty quantification (UQ) setting. Building directly on the neural cDFT formulation, we interpret the equilibrium one-body density field as a physical likelihood, exploiting its statistical-mechanical meaning as the probability density of particle positions. Molecular simulation trajectories are incorporated directly as observations, and inference is performed over the parameters governing the learned free-energy corrections. A key enabling component of the proposed methodology is the adjoint formulation for nonlocal operators developed in the underlying neural cDFT framework. This adjoint construction provides efficient gradients of the equilibrium density with respect to free-energy parameters, despite the implicit and nonlocal structure of the cDFT Euler--Lagrange equations. As a result, it enables the practical application of modern Bayesian inference techniques---including variational inference with normalizing flows and Hamiltonian Monte Carlo---while treating the equilibrium density field itself as the likelihood function. The resulting Bayesian formulation enforces the cDFT Euler--Lagrange equations as hard constraints, ensuring that all posterior samples correspond to admissible thermodynamic states. The framework yields full posterior distributions over free-energy corrections, equilibrium density fields, and derived macroscopic observables, enabling principled uncertainty propagation from molecular data to continuum-scale predictions. This work establishes a scalable route to uncertainty-aware surrogates for molecular fluids predictions, combining statistical mechanics, adjoint method, machine learning, and modern UQ.