Edge-on protoplanetary disks provide a spatially resolved view of the vertical distribution of icy planetary building blocks, but inferring their relative abundances requires radiative-transfer (RT) modeling with multiple anisotropic scattering. Each model costs 10¹-10³ CPU hours, which precludes principled Bayesian inference and homogeneous, population-scale analyses. We develop and validate a fast and reusable surrogate for Bayesian inference with ice-bearing edge-on disk models. We train a Fourier Neural Operator on a grid of RADMC-3D spectra. The network maps disk parameters and analytic physical input fields to the spectrum, using mixed-resolution spectra. We use the surrogate in an explicit-likelihood nested sampler and validate it on held-out simulations and synthetic recovery tests with TARP coverage. Finally, we use it to retrieve the disk structure and ice abundances of ESO-Hα 569 from its spectral energy distribution. The surrogate reproduces the RADMC-3D spectra with a median error of 4.84%. We recover synthetic parameters without detectable bias and with nominal coverage. For ESO-Hα 569, the posterior favors a massive (M_(rm gas)=0.028^(+0.015)_(-0.011) M_⊙), nearly edge-on (i=79.8^(+5.2)_(-2.4) deg) disk consistent with the image-plus-SED model of Wolff et al. (2017). Its ice abundances and cavity size remain weakly constrained by broadband data alone. The marginal spectrum-evaluation speed-up is ∼10⁹-10¹⁰. Neural-operator emulation makes full posterior inference feasible within the adopted model and training prior. Key posterior points should still be checked with high-fidelity RT models. By drastically cutting evaluation time, this approach enables for the first time homogeneous population studies of ices in disks, which new JWST surveys can capitalize on to map the range of compositions within planetary nurseries.