Variability is a universal feature among biological units such as neuronal cells as they enable a robust encoding of a high volume of information in neuronal circuits and prevent hypersynchronizations such as epileptic seizures. While most computational studies on electrophysiological variability in neuronal circuits were done with simplified neuron models, we focus on the variability of detailed biophysical models of neurons. Using measures of experimental variability, we leverage a Markov chain Monte Carlo method to generate populations of electrical models that can reproduce the variability from sets of experimental recordings. These models are constrained via a cost function to match specific electrical features extracted from experimental recordings. We demonstrate our approach on layer 5 pyramidal cells with continuous adapting firing type and show, for example, that morphological variability is insufficient to reproduce electrical variability. Not only does this method enable statistical analysis of neuronal models, but it also helps build accurate models, by detecting consistent biases from the data. Overall, this approach provides a strong statistical basis to create detailed models of neurons with controlled variability to study specific dendritic mechanisms such as nonlinear dendritic computations in the apical dendrites, bursting in thalamic cells and eventually help understand the impact of single-cell variabilities on neural circuit dynamics, such as synchronisation properties.