The flux variability of blazars at very high energies does not have a clear origin. Flux variations on time scales down to the minute suggest that variability originates in the jet, where a relativistic boost can shorten the observed time scale, while the linear relation between the flux and its rms or the skewness of the flux distribution suggests that the variability stems from multiplicative processes, which are associated in some models with the accretion disk. We study the rms-flux relation and emphasize its link to Pareto distributions, characterized by a power-law probability density function. Such distributions are naturally generated within a minijets-in-a-jet statistical model, in which boosted emitting regions are isotropically oriented within the bulk relativistic flow of a jet. We prove that, within this model, the flux of a single minijet is proportional to its rms. This relation still holds when considering a large number of emitting regions, for which the distribution of the total flux is skewed and could be interpreted as being log-normal. The minijets-in-a-jet statistical model reconciles the fast variations and the statistical properties of the flux of blazars at very high energies.