  Type-I X-ray bursts are recurring thermonuclear explosions on the surface of
accreting neutron stars. Matching observed bursts to computational models can
help to constrain system properties, such as the neutron star mass and radius,
crustal heating rates, and the accreted fuel composition, but systematic
parameter studies to date have been limited. We apply Markov chain Monte Carlo
methods to 1D burst models for the first time, and obtain system parameter
estimations for the `Clocked Burster', GS 1826$-$238, by fitting multiple
observed epochs simultaneously. We explore multiple parameters which are often
held constant, including the neutron star mass, crustal heating rate, and
hydrogen composition. To improve the computational efficiency, we precompute a
grid of 3840 KEPLER models - the largest set of 1D burst simulations to date -
and by interpolating over the model grid, we can rapidly sample burst
predictions. We obtain estimates for a CNO metallicity of $Z_\mathrm{CNO} =
0.010^{+0.005}_{-0.004}$, a hydrogen fraction of $X_0 = 0.74^{+0.02}_{-0.03}$,
a distance of $d \sqrt{\xi_\mathrm{b}} = 6.5^{+0.4}_{-0.6}\, \mathrm{kpc}$, and
a system inclination of $i = {69^{+2}_{-3}}^{\circ}$.
