We develop a machine-checked construction of non-Gaussian functional integrals using the Henstock–Kurzweil gauge integral and Chernoff product approximations. The central object is a finite family of bosonic modes with action S(phi) = (1/2) phi^T A phi + lambda * sum_i phi_i^4, where A is positive definite. We prove that the one-mode integral I(omega, j, lambda) is finite, strictly positive, monotone and infinitely differentiable in the coupling lambda on [0, infinity). Its derivatives are given by convergent integrals of phi^{4k} with the same weight, not by the divergent perturbative series. The M-mode influence functional factorises into one-mode integrals and is bounded by its Gaussian value. A Chernoff / Lie–Trotter splitting handles the non-commutativity of the free and non-Gaussian generators. All statements are formalised in Lean 4 with Mathlib; the accompanying file HkNonGaussian.lean is free of sorry and uses only the standard axioms propext, Classical.choice, Quot.sound. Four illustrative applications are worked out at the level of explicit formulas: the Duffing oscillator, local volatility (CEV) in finance, Wilson–Cowan neural fields, and non-Gaussian quantum reservoirs. The construction is completely direct and does not use Wick rotation, Wiener measure, zeta-regularisation or analytic continuation back from imaginary time.