Random number generation from a power law

General Usage
,
1

hello all,
I am new to Julia. I need to generate random numbers from power law, broken power law and different distributions. Is there a package that do all of these? Thank you very much in advance.

2

Welcome to the discussion board. :slight_smile:
I have not used this, but perhaps Distributions.jl might work for you?

3

More precisely, here https://juliastats.org/Distributions.jl/stable/univariate/#Sampling-(Random-number-generation) you can find the univariate distributions and sampling from an arbitrary one.
Example:

using Distributions, Plots
plotlyjs()
const α = 2.75
const xmin= 3.0
d = Pareto(α, xmin)
smpl= rand(d, 1000)
h=histogram(smpl, bins= xmin:50, normalize=:pdf, 
            framestyle=:box, size=(500, 250), legend=false)

histPareto
histPareto700×500 6.08 KB

Broken power law distribution is not implemented in Distributions.jl. It depends on the number of break points and power indices on each interval.
I only know of the existence of the broken power law, implemented in astropy: https://docs.astropy.org/en/stable/_modules/astropy/modeling/powerlaws.html, but without a sampling method.

4

Don’t know much about broken power laws, but mixture models might be able to piece them together:

alpha_below = 1.6
alpha_above = 1.2
breakpoint = 5.4
p_below = cdf(Pareto(alpha_below), breakpoint)
broken_dist = MixtureModel([truncated(Pareto(alpha_below), nothing, breakpoint),
                            truncated(Pareto(alpha_above), breakpoint, nothing)],
                           [p_below, 1 - p_below])
rand(broken_dist)
1 Like
5

If your distribution is known and sufficiently smooth[1] you can use ApproxFun.jl.

  1. Though in practice 2M Chebyshev modes is enough to approximate some very wonky functions. ↩︎

1 Like
6

I implemented broken-power-law sampling in InitialMassFunctions.jl, which is now a registered package. This functionality is provided by the BrokenPowerLaw type. This type implements the Distributions.jl sampler API, so a more efficient object for sampling can be constructed via Distributions.sampler(::BrokenPowerLaw). You may also be interested in the single power-law PowerLawIMF type, which is just a wrapper on top of the Pareto distribution from Distributions.jl.

1 Like
7

All of this is really helpful. Thanks a lot!