Extreme TeV BL Lacs are a class of blazars with unique spectral and temporal features that are not easily reproducible using standard one-zone models based on single shock acceleration. To account for their peculiar properties, we elaborated a two-step acceleration model in which a recollimation shock and the subsequent downstream turbulence energize non-thermal electrons. We applied the model to a sample of extreme TeV BL Lacs with well-characterized spectral energy distributions. Since we used several sources, we automatized the exploration of the parameter space. This allowed us to derive the parameter distributions and study the correlations among them. We numerically solved a system of two coupled nonlinear differential equations to obtain the non-thermal particles and turbulence spectra. We calculated the spectral energy distribution via the synchrotron self-Compton emission model. The automatization of the parameter space exploration is possible through a Markov chain Monte Carlo (MCMC) ensemble sampler, in our case emcee We derived well-defined posterior distributions for the parameters, showing that the model is well constrained by available data and demonstrating the suitability of our method. The cross-correlations among some of the physical parameters are not trivial.\ Therefore, we conclude that MCMC sampling is a key instrument for characterizing the complexity of our multiparameter phenomenological model.