This study introduces Bayesian neural networks (BNNs) as a model of nonlinear flame dynamics. Other than standard neural networks, BNNs can quantify uncertainty in a principled manner. Such uncertainty estimation is crucial for the reliable modeling of thermoacoustic combustion instabilities. To determine an ideal network architecture, we employ evidence-based model selection. This approach minimizes model size and complexity, thereby reducing epistemic uncertainty without compromising prediction quality.\\ The framework is applied to a well-studied laminar premixed methane-air slit flame. BNNs are trained on broadband-forced incompressible CFD simulation data and tested on both broadband and harmonically forced data. Analysis of the resulting uncertainties demonstrates that BNNs can effectively indicate a lack of data or physical information. The predicted flame describing function uncertainties (Figure \ref{fig:FDF_06s}) confirm previous reports that broadband-forced flames respond less nonlinearly than harmonically forced flames with similar forcing amplitudes. Furthermore, uncertainty quantification indicates that time series data generated with nonlinear broadband excitation contains only little linear information. Finally, the optimized BNN is coupled with an acoustic solver to accurately predict the limit cycle of an intrinsic thermoacoustic instability.