We present a method to self-consistently propagate stellar-mass [M⋆=log(M/M⊙)] and star-formation-rate [Ψ=log(ψ/M⊙yr−1⁠)] uncertainties on to intercept (α), slope (β), and intrinsic-scatter (σ) estimates for a simple model of the main sequence of star-forming galaxies, where Ψ=α+βM⋆+N(0,σ)⁠. To test this method and compare it with other published methods, we construct mock photometric samples of galaxies at z ∼ 5 based on idealized models combined with broad- and medium-band filters at wavelengths 0.8–5 μm. Adopting simple Ψ estimates based on dust-corrected ultraviolet luminosity can underestimate σ. We find that broad-band fluxes alone cannot constrain the contribution from emission lines, implying that strong priors on the emission-line contribution are required if no medium-band constraints are available. Therefore, at high redshifts, where emission lines contribute a higher fraction of the broad-band flux, photometric fitting is sensitive to Ψ variations on short (∼10 Myr) time-scales. Priors on age imposed with a constant (or rising) star formation history (SFH) do not allow one to investigate a possible dependence of σ on M⋆ at high redshifts. Delayed exponential SFHs have less constrained priors, but do not account for Ψ variations on short time-scales, a problem if σ increases due to stochasticity of star formation. A simple SFH with current star formation decoupled from the previous history is appropriate. We show that, for simple exposure-time calculations assuming point sources, with low levels of dust, we should be able to obtain unbiased estimates of the main sequence down to log(M/M⊙)∼8 at z ∼ 5 with the James Webb Space Telescope while allowing for stochasticity of star formation.