  Power-law distributions are essential in computational and statistical
investigations of extreme events and complex systems. The usual technique to
generate power-law distributed data is to first infer the scale exponent
$\alpha$ using the observed data of interest and then sample from the
associated distribution. This approach has important limitations because it
relies on a fixed $\alpha$ (e.g., it has limited applicability in testing the
{\it family} of power-law distributions) and on the hypothesis of independent
observations (e.g., it ignores temporal correlations and other constraints
typically present in complex systems data). Here we propose a constrained
surrogate method that overcomes these limitations by choosing uniformly at
random from a set of sequences exactly as likely to be observed under a
discrete power-law as the original sequence (i.e., regardless of $\alpha$) and
by showing how additional constraints can be imposed in the sequence (e.g., the
Markov transition probability between states). This non-parametric approach
involves redistributing observed prime factors to randomize values in
accordance with a power-law model but without restricting ourselves to
independent observations or to a particular $\alpha$. We test our results in
simulated and real data, ranging from the intensity of earthquakes to the
number of fatalities in disasters.
